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SUMMARY 

This  report  describes  the  results  of  our  efforts  on  AFOSR  Contract  F49620- 
77-C-0076,  Numerical  Model  Development  for  Laser  Cavity  Flowfields.  The  research 
was  aimed  toward  developing  a  technique  for  automatically  varying  (in  time)  the 
one-dimensional  spatial  finite  difference  mesh  over  which  certain  time  dependent 
partial  differential  equations  were  discretized.  In  particular,  the  variable 
mesh  points  should  adapt  themselves  to  the  solution  that  is  being  marched  out  in 
time,  i.e.,  they  should  be  closely  spaced  in  regions  where  there  are  large 
spatial  gradients  and  widely  spaced  where  the  solution  is  smooth.  Further,  this 
automatic  mesh  variation  method  is  to  be  designed  to  be  used  in  conjunction  with 
the  GEARIB  implicit  integration  package,  which  was  used  successfully  for  a  variety 
of  fixed  mesh  discretizations  during  the  first  two  years  of  this  contract. 

(GEARIB  is  used  march  out  in  time  the  solution  of  the  system  of  ordinary 
differential  equa  ons  resulting  from  spatially  discretizing  the  governing  PDE(s).) 

The  partial  differential  equation  that  was  used  as  a  test  bed  for  this 
automatic  mesh  variation  was  Burger's  equation.  Numerical  techniques  (using 
GEARIB  on  a  fixed  mesh)  for  solving  this  PDE  were  developed  during  the  first  year 
of  this  contract.  Although  Burger's  equation  is  simple  in  appearance,  its 
solutions  have  many  of  the  characteristics  (nonl inearity,  wave-like  propagation, 
diffusive  effects,  and  steep  gradients)  that  are  found  in  solutions  of  laser  flow 
equations,  which  are  the  problems  that  the  methods  developed  herein  are  ultimately 
intended  for. 
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I.  ADAPTIVE  MESH  FOR  BURGER'S  EQUATION 


A.  Introduction 

In  the  report  on  the  results  of  the  first  year  of  this  study.  Reference  [1], 
a  detailed  discussion  was  given  in  Section  II  regarding  the  choice  of  an 
integration  method  to  march  out  the  solution  in  the  time-like  coordinate.  It  was 
concluded  that  an  implicit  integration  method  would  be  required  to  efficiently 
handle  differential  systems  containing  the  strong  stiffness  effects  due  to 
chemistry  and  radiation  present  in  the  numerical  solution  of  laser  flows .  It  was 
also  concluded  that  the  GEARIB  package  [2],  which  embodies  a  family  of  implicit 
methods  and  which  is  also  designed  to  treat  implicit  differential -algebraic 
systems  was  well  suited  to  our  applications. 

During  the  second  year  of  this  study,  numerical  difficulties  arose  in  the 
treatment  of  the  laser  flow  equations  because  the  onset  of  lasing  was  accompanied 
by  extreme  spatial  gradients  in  some  of  the  dependent  variables  being  integrated. 

The  location  of  these  gradients  is  unknown  a  priori,  and  may  in  general  vary  in 
time.  In  order  to  obtain  solutions,  a  trial  and  error  mesh  selection  process  was 
used;  pick  a  mesh,  solve  the  problem,  refine  the  mesh  based  on  the  current  solution, 
solve  the  problem  again,  etc.  The  only  way  to  avert  such  an  iterative  procedure 
would  be  to  use  a  uniformly  fine  spatial  mesh.  Such  a  mesh  adequate  to  accommodate 
the  gradients  of  realistic  laser  flows  would  be  computationally  prohibitive.  Thus, 
it  became  apparent  that  a  successful  numerical  treatment  of  such  flows  would 
require  a  finite  difference  mesh  that  automatically  adapts  to  the  solution. 

This  requirement  for  an  adaptive  spatial  finite  difference  mesh  has  been 
cited  previously  in  the  literature.  For  example,  Dwyer  and  Sanders  [4]  in  their 
study  of  unsteady  combustion  phenomena  suggest  an  approach  to  an  adaptive  finite 
difference  mesh  but  do  not  actually  use  it.  (All  the  results  they  present  are 
for  uniform  meshes.)  Chong  [5]  uses  adaptive  meshes  in  his  study  of  Burger's 
equation,  which  has  solutions  with  steep  spatial  gradients.  Although  he  obtains 
good  results,  his  approach  does  not  seem  practical  for  our  applications.  The 
main  reason  is  that  Chong  adjusts  his  finite  difference  mesh  at  discrete  times, 
which  would  mean  in  our  case  that  the  GEARIB  integrator  would  have  to  be  restarted 
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after  each  such  adjustment.  This  is  a  rather  severe  penalty  to  pay  especially  if 
GEARIB  had  been  operating  with  a  high  order  method  prior  to  the  rediscretization, 
since  restart  begins  with  a  first  order  method  and,  consequently,  a  rather  small 
step  size.  A  similar  approach  has  been  taken  by  Miller,  Morton,  and  Baines  [6]  in 
a  moving  boundary  problem  governed  by  a  generalized  one-dimensional  diffusion 
equation. 

In  the  approach  described  in  this  report,  the  rediscretization  is  done  on  a 
continuous  basis,  and  so  there  is  no  need  to  restart  the  integrator  (as  there  is 
when  the  spatial  mesh  is  altered  at  discrete  times).  A  further  advantage  of  this 
approach  is  that  moving  boundary  problems  are  readily  accommodated. 

B.  Discretization  of  Burger's  Equation 

One  of  the  simplest  nonlinear  PDE's  which  contains  both  convective  and 
diffusive  effects  is  Burger's  equation. 


IT 

at 


32T 

T? 


t  9T 


(1) 


where  B  and  y  are  constants  which  can  be  varied  to  emphasize  the  diffusive  or 
convective  contributions.  Because  Burger's  equation  has  known  exact  solutions 
which  are  traveling  waves  with  very  steep  fronts  (like  shock  waves)  it  has  been 
a  popular  equation  to  test  numerical  schemes  on.  In  fact,  during  the  first  year 
of  this  study  a  detailed  numerical  treatment  of  this  equation  was  undertaken 
(see  Chapter  IV  of  Schimke,  Rushmore,  and  Zelazny  [1]  for  details).  Because  of 
this  experience  with  Burger's  equation,  it  seems  an  ideal  test  bed  for  our 
proposed  adaptive  mesh  scheme. 

As  in  Reference  [1],  we  begin  by  introducing  the  auxiliary  variable,  v. 


v 


_  9T 
"  3x 


Then  (1)  can  be  replaced  by 


3T  _ 
3t 


B  f  -  yvT 


(2) 


(3a) 
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The  boundary  conditions  are  assumed  to  be  separated  but  otherwise  quite 
general ,  i.e., 


9-|  (^1  ,V1  *t)  ®N^N*VN’^  ®  (4) 

where  g]  and  g^  are  arbitrary  functions  continuous  in  t,  and  the  subscripts  1 
and  N  refer  to  the  left  and  right  boundaries,  which  may  be  moving  in  time. 

Now  assume  a  finite  difference  mesh  with  nodes  at  x-j ,  x2»  ...»  x^  (which 

are  in  general  functions  of  time),  with  and  associated  with  x^,  and 

H.e(x.+i  -  x.j ) .  Now  since  T..  corresponds  to  a  mesh  point  which  is  moving,  the 

rate  of  change  of  T.  in  time  (i.e.,  T . )  is  not  the  same  as  [3T/3tl  in  (3a). 

11  x. 

The  relation  between  these  two  quantities  is: 


Now  the  x.j  are  to  be  determined  automatically  during  the  course  of  the 
problem  so  we  must  establish  conditions  which  will  uniquely  specify  them.  For 
x-|  and  x^  we  assume. 


X1  =  xn  =  V1^ 

where  f-j  and  f^  are  prescribed  functions  of  time,  f-j  and  fN  could  be  readily 
generalized  to  be  functions  of  the  dependent  variables  without  complicating  the 
numerical  process.  Such  f^  and  f^  are  required  in  certain  types  of  moving 
boundary  problems,  but  we  will  not  consider  them  here. 

To  complete  the  system  of  equations,  relationships  specifying  the  interior 
x.j  must  be  established.  This  specification  is  the  real  crux  of  the  numerical 
problem.  We  certainly  expect  that  a  suitable  mesh  spacing  should  vary  inversely 
with  the  local  gradients,  i.e.,  the  local  values  of  v^  Such  a  relationship  can 
be  written: 


(xi+1  -  x.)  = 


F(vi’vi+1} 


(7) 
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where  Fisa  function  that  can  be  specified  in  various  ways  and  A  is  a  positive 
function  of  t  but  is  constant  in  x.  A  in  effect  provides  a  uniform  scaling  of 
all  the  (x.j+i  -  x.j),  i  =  l.N-l  so  that  the  mesh  "fits  into"  the  region  [x-j.x^]. 

F  in  (7)  could  have  been  made  a  function  of  T\  and  T.+-j  (in  addition  to  v.  and 
vi+-|)  without  introducing  any  essential  complication,  but  we  have  had  no  occasion 
to  do  so  during  the  course  of  this  study.  The  relation  (7)  involves  values  at 
two  adjacent  mesh  points.  Why  not  allow  a  more  general  relation  involving  three 
(or  more)  consecutive  mesh  points?  This  would  be  contrary  to  the  philosophy  which 
led  to  the  introduction  of  the  auxiliary  dependent  variable,  v,  in  equation  (2). 
Introduction  of  v  allows  a  simple  differencing  scheme  to  be  employed  for  the 
partial  differential  system  (10).  This  differencing  scheme  is  unaffected  by  mesh 
nonuniformity,  requires  no  modification  next  to  boundary  points,  and  results  in 
a  Jacobian  matrix* with  a  narrow  bandwidth. 

According  to  equation  (7),  (x^  -  x^)  is  given  explicitly  as  a  function  of 
A,  v^  and  vi+1.  In  practice,  it  is  more  convenient  to  specify  (xi+1  - 
implicitly.  Instead  of  (7)  we  use 


{xi 


+1 


*i> 


F(xi’  xi+1,  vi ,  v.+1) 


A  specific  example  of  (8)  used  in  the  course  of  this  study  is: 


(8) 


A 


2  2 

The  squared  term  in  (9)  is  an  approximation  to  •  T/<x  .  Thus,  (9)  requires  that 

2  2 

the  mesh  become  finer  in  regions  of  large  "■  T/  >x  ,  a  reasonable  criterion.  The 

positive  constant  e  is  present  to  prevent  the  mesh  spacing  from  becoming  arbitrarily 

2  2 

large  in  regions  of  small  ,»  T/>x  .  Note  that  F(xi ,  xi+-j,  vi ,  vi+-j)  as  defined  in 
(9)  is  always  positive,  thereby  guaranteeing  that  consecutive  mesh  points  will 
always  be  distinct.  (A  can  never  vanish.) 


* I mp licit  integration  methods  (as  in  GEARIB)  require  the  formation  and  factorization 
of  the  Jacobian  matrix  of  the  differential  system. 
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When  A  was  introduced  in  (7),  it  was  described  simply  as  a  scaling  factor. 

But  it  also  has  a  "geometric"  interpretation.  Rewriting  (8)  and  summing  over 
the  (N-l)  mesh  intervals  between  x-j  and  we  get: 

i  N_1  i  f*N 

A  "  IfTTT  i^]  F(xi’xi+l’vi,vi+Vxi+l  '  V  ~  j  F^xi  )Xi+i  >vi  >vi+l ^dx  (10) 

X1 

Now  F(x .  >x.j+i  ,v^  >v -+i )  is  normally  selected  to  be  some  measure  of  the  "variability" 
of  the  solution  within  [x^x^].  For  example,  in  (9)  as  e->0 

I  2  I*5 

r  -  U  T  i 


Thus,  A  is  a  measure  of  the  average  variability  over  a  mesh  interval. 

Now  differencing  the  partial  differential  system  (3)  in  an  obvious  way, 
using  (4),  (5),  (6),  and  (8),  and  adopting  the  notation. 


Hi  xi+l  -  *1 

(12a) 

7io  5  (,i  +  v0)/z 

(12b) 

Tu  5  <Ti*y/2 

(12c) 

we  obtain  the  differential -al gebraic  system  on  the  following  page. 

Notice  that  product  terms  such  as  vT  in  (3a)  and  xv  in  (5)  have  been 
approximated  as  products  of  sums  in  (13)  (instead  of  sums  of  products).  For 
example,  on  mesh  interval  i,  we  wrote  the  term  vT  as 

7i,i+i  Ti,i+i  =  (vi  +  vi+i)(Ti  +  W/4  instead  of  (Vi  +  Vi  W/2- 
Numerical  experiments  were  performed  using  both  formulations.  The  results, 

although  not  conclusive,  tended  to  favor  the  product  of  sums  formulation. 

Possibly  the  reason  is  that  the  effect  of  spurious  high  frequency  oscillations 

(especially  in  v.)  tend  to  be  negated  when  using  products  of  sums. 


1-5 


Bell  Aerospace 


TEXTRON 


Although  the  matrix  multiplying  the  derivative  vector  is  singular,  this 
presents  no  problems  for  the  GEARIB  integrator.  (Similar  singular  differential  - 
algebraic  systems  were  solved  during  the  first  two  years  of  this  study.)  One 
special  characteristic  of  system  (13)  which  was  not  present  in  the  numerical 
treatment  of  Burger's  equation  on  a  fixed  mesh,  is  that  the  Jacobian  of  the  system 
is  no  longer  strictly  banded.  The  last  column  has  many  non-zero  elements.  Also, 
certain  variants  of  system  (13)  having  many  non-zero  elements  in  the  last  row  of 
the  Jacobian  have  been  tested.  Consequently,  it  was  necessary  to  modify  GEARIB 
(which  is  tailored  to  systems  with  strictly  banded  Jacobians)  to  account  for 
the  possibly  full  last  row  and  column.  The  GEARIB  LU  decomposition  routine  was 
modified  so  that  the  upper  (N-l)x(N-l)  is  treated  as  a  band  matrix  and  special 
formulas  are  used  for  the  last  row  and  column  of  the  L  and  U  factors,  respectively. 
With  this  modified  LU  decomposition  plus  a  few  changes  within  some  of  the  other 
routines  of  the  GEARIB  package,  such  systems  with  "almost"  banded  Jacobians  are 
efficiently  integrated. 

C.  Specification  of  Boundary  and  Initial  Conditions 

In  Section  IV.C.l  of  [1]  an  extensive  discussion  is  given  regarding  the 
specification  of  boundary  and  initial  conditions  when  solving  Burger’s  equation 
on  a  fixed  mesh.  It  is  shown  that  it  ’  s  important  to  choose  initial  and  boundary 
conditions  so  that  all  strictly  algebraic  equations  are  satisfied  initially.  On 
a  fixed  mesh,  this  is  easy  to  do  for  the  boundary  conditions 


=  g^t). 

Tn  =  gN(t),  or 

(14a) 

=  g^t). 

VN  ■  9N(t) 

(14b) 

Then  the  algebraic  equations  arising  from  the  discretization  of  (3b),  i.e., 

■  <Ti  *  W/Hi  (’5) 

are  satisfied  by  solving  (15)  recursively  for  vi  (assuming  T.  is  prescribed  which 
is  the  usual  case).  For  boundary  conditions  other  than  (14),  the  initial 
satisfaction  of  all  algebraic  equations  is  less  straightforward.  The  additional 
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complications  arising  out  of  such  boundary  conditions  and  methods  for  dealing 
with  them  are  detailed  in  [1]. 

For  variable  meshes,  there  is  the  additional  set  of  algebraic  equations 
obtained  by  writing  (8)  at  the  mid  points  of  each  mesh  interval,  i.e.,  equations 
(5),  (8),  ...,  (N-2)  of  system  (13).  These  equations  (e.g.,  (9))  are  nonlinear. 
Consequently,  their  initial  solution  poses  more  of  a  problem.  In  view  of  this, 
all  of  the  examples  described  in  this  report  use  boundary  conditions  (14a)*  in 
order  not  to  compound  the  problem  of  determining  initial  conditions  and  therefore 
deter  us  from  our  main  goal  of  devising  a  reliable  scheme  for  automatically 
adapting  the  finite  difference  mesh  to  the  solution  as  it  evolves  in  time. 

Another  assumption  that  is  made  in  all  the  examples  of  this  report  is  that 
the  initial  values  of  T  are  given  by  a  function  of  x  which  is  continuous  and 
piecewise  differentiable.  The  method  for  determining  initial  conditions  for 
x.j  and  v.j  described  in  the  following  is  based  on  this  assumption.  (A  discussion 
of  the  problem  where  the  initial  T,.  are  assigned  values  at  discrete  x^  is  given 
at  the  end  of  this  section.) 

A  modified  Newton  method  was  used  to  solve  the  system  of  nonlinear  equations 
associated  with  the  determination  of  the  initial  conditions.  The  GEARIB  package 
is  designed  to  solve  differential -algebraic  systems,  which  includes  purely 
algebraic  systems.  Thus,  the  solution  for  the  initial  conditions  can  be  done 
within  the  framework  of  GEARIB.  Certain  changes  to  the  standard  GEARIB  scheme 
were  found  useful  to  increase  the  likelihood  of  convergence. 

For  the  determination  of  the  initial  conditions,  the  governing  system  is  (13) 
with  the  matrix  on  the  LHS  set  to  zero  and  with  equations  (3),  (6),  ...,  (N-4) 
replaced  by 


0 


Ti 


T(x.) 


06) 


*(14a)  was  chosen  over  (14b)  since  our  sample  problems  all  have  wave-like  solutions 
which  travel  in  the  positive  x-direction.  Consequently,  the  "upstream"  boundary 
condition  tends  to  have  a  less  important  effect  on  the  solution,  especially  when  the 
wave  impinges  on  the  downstream  boundary.  Thus,  we  have  prescribed  the  more  usual 
condition  (on  T  rather  than  v)  at  the  downstream  boundary. 
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where  T(x)  is  the  function  defining  the  initial  values  of  T  as  a  continuous 
function  of  x.  To  distinguish  the  iterations  for  the  initial  conditions  within 
GEARIB,  an  artificial  time  starting  with  -50  and  incremented  by  1  on  each  call 
to  GEARIB  was  established.  (Integration  of  the  actual  differential  system 
began  at  t=0.)  This  allows  50  iterations  for  the  initial  conditions  to  converge. 
Convergence  typically  took  place  within  3  to  12  iterations.  The  problem  was 
aborted  if  50  iterations  did  not  produce  convergence. 


GEARIB  requires  a  number  of  user  coded  routines  to  define  the  system  to  be 
solved.  Within  each  of  these  routines  there  is  at  least  one  test  on  t.  If  t 
is  negative,  branches  are  taken  to  define  the  system  for  the  initial  conditions; 
for  t  positive,  branches  are  taken  to  define  system  (13).  Error  tests,  convergence 
tests,  and  the  Jacobian  re-evaluation  criterion  in  GEARIB  are  modified  so  that 
only  one  Newton  Iteration  (with  a  newly  evaluated  Jacobian)  is  performed  for  each 
call  to  GEARIB  when  solving  for  initial  conditions.  (Testing  for  convergence 
of  the  iterations  is  performed  in  the  calling  program.) 


The  specific  steps  of  the  initial  conditions  iteration  procedure  are  as 
follows: 


1.  Set  time  to  -50  and  step  size  to  1. 

2.  Choose  x.  to  be  equally  spaced  on  the  interval  [x^x^]. 

3.  Evaluate  T-  from  (16). 

4.  Solve  for  v^  recursively  using  (15). 

5.  Solve  for  A  using 


where  F  is  as  in  (8).  (Alternatively,  equation  (10)  could  be  used  for  A.) 

6.  Place  the  values  of  all  the  dependent  variables  (T..,  v^ ,  xi ,  A)  in  a  "save" 
vector. 


7.  Call  GEARIB  to  do  one  Newton  iteration. 
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8. 


9. 


10. 


11. 

12. 


Check  to  see  whether  any  mesh  interval  (xi+-j  -  x^)  has  decreased  by  more 
than  20%  during  the  Newton  iteration.  If  so,  relax  the  Newton  iteration 
by  performing  a  uniform  linear  interpolation  (for  all  x^)  between  the 
current  and  preceding  values  so  that  the  maximum  decrease  in  any  mesh 
interval  is  20%. 

Evaluate  T.  from  (16). 

Solve  for  v^  recursively  using  (15). 

Solve  for  A  as  in  Step  (5). 


Test  for  convergence,  i.e.. 


3N+1 


[(3N+1)EPS]2 


? 


where 


i)  The  summation  is  over  all  the  dependent  variables  (Ti#  vi ,  xi ,  A). 

ii)  AY.  is  the  difference  between  the  current  and  preceding  (at  the 

J  t  h 

beginning  of  Step  (7))  values  of  the  j  dependent  variable. 

iii)  YN.  is  the  "normalization"  factor  for  dependent  variable  j.  It 

J 

is  the  maximum  value  of  Y •  occurring  during  all  of  the  preceding 
iterations,  but  it  must  be  at  least  0.01. 

iv)  EPS  is  the  same  accuracy  criterion  used  for  the  conventional  GEARIB 
integration. 


13.  If  the  iteration  has  converged  go  to  the  beginning  of  the  conventional 
integration. 


14.  Test  to  see  if  the  maximum  number  of  iterations  (=50)  has  been  exceeded 
or  if  there  is  an  error  return  from  GEARIB.  If  so,  then  abort. 


15.  Go  to  Step  (6) . 


As  mentioned  previously,  the  principal  focus  of  the  work  described  in  this 
report  was  the  study  of  how  the  variable  mesh  adapts  itself  to  the  evolving  solution. 
Consequently,  the  emphasis  in  the  initial  condition  iteration  scheme  was  on  the 
assurance  of  convergence  rather  than  on  efficiency.  This  is  the  reason  for  the 
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conservative  relaxation  criterion  of  Step  8  and  the  re-evaluation  of  the  Newton 
matrix  on  each  iteration.  Despite  these  safety  factors,  the  time  required  for  the 
determination  of  initial  conditions  was  generally  an  insignificant  fraction  of  the 
total  run  time. 

What  about  the  commonly  occurring  situation  where  the  initial  values  of  T\ 
are  prescribed  at  preselected  values  of  x^?  The  algebraic  equations  (15)  could 
be  satisfied  as  before  by  solving  for  the  v^  recursively.  However,  equations  (8) 
could  not  be  satisfied.  A  way  out  of  this  difficulty  would  be  to  generalize  (8) 
so  that  the  function  F  would  be  defined  somewhat  differently  on  each  mesh  interval, 
but  in  such  a  manner  that  -*  F  with  increasing  time.  If  F  specifies  the 

"optimal"  distribution  of  the  xi ,  then  initially  and  for  some  time  the  x^ 
distribution  actually  used  would  be  suboptimal ,  but  this  cannot  be  avoided. 

All  of  the  F  definitions  used  during  this  study  contain  a  parameter  e,  e.g., 
equation  (9).  An  easy  way  to  implement  the  scheme  just  described  would  be  to 
preselect  A  and  then  replace  e  by  c.Q  so  that  (8)  would  be  exactly  satisfied  on 
each  interval.  Then  define 

ei  =  ei()e"at  +  e(l  -  e'at) 

where  a  is  some  positive  constant  to  be  determined;  probably  by  experiment. 

D.  Specification  of  Finite  Difference  Mesh  Spacing 

As  noted  previously,  the  principal  difficulty  of  this  study  is  the  selection 
of  the  function  F(x . ,x . . , ,v . ,v )  of  equation  (8).  To  motivate  this  selection 

1  I  T  I  1  1  +  1 

consider  the  spatial  discretization  error  when  (3a)  is  written  at  the  mid  point 
of  each  mesh  interval,  i.e.,  equations  (3),  (6),  ...,  (N-4)  of  system  (13). 

‘A  *  W  -  <vi  *  +  W/2 

■  26<Vl  -  vp/H,  -  Y(V,  *  »1+1)(T,  ♦  Tj+1)/2  (17) 

Now  expand  each  term  of  (17)  about  the  midpoint  cf  the  mesh  interval.  Quantities 
associated  with  this  midpoint  are  given  a  subscript  m.  Apostrophes  denote  9/ ax. 
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.  h*  . 
2t  +  r-  r 

m  4  m 


H?  \/.  H?  . 

*  r  +  r  xro 
2ek  +  2T  vm")  ' 


■<(2\ +  r  v)(Tm  *  r  Tm)  +  °(Ht ) 


Thus,  the  truncation  error  associated  with  this  spatial  discretization  is: 


4  (3  vn 


-  y(v  T"  +  T  v")  -  T"  +  (v  x"  +  x  v'’")/  +  O(H^)  (19) 

''mm  m  m  '  m  'mm  mmy/  '  V  '  ' 


Note  that  the  truncation  error  is  0(hr)  irrespective  of  any  mesh  nonuniformity. 

Contrast  this  to  the  classical  three  point  differencing  scheme  for  which  the 

2 

truncation  error  is  0(H  )  only  for  uniform  meshes.  Thus,  for  a  preselected 
truncation  error  level  TE,  should  be  chosen, 


Hi  = 


2/fE 

/rrn 


There  are  several  reasons  why  (20)  is  impractical: 

1)  To  establish  a  finite  difference  mesh  satisfying  (20),  it  is  necessary 

to  fix  at  least  one  mesh  point,  say  x-j .  Then  assuming  that  {  }  can 

be  evaluated,  (20)  can,  in  principle,  be  used  to  specify  all  the  mesh 
points.  However,  for  an  arbitrarily  selected  value  of  TE,  none  of  the 
mesh  points  so  determined  would,  in  general,  coincide  with  the  desired 
right  boundary  point,  x^.  Even  if  TE  could  be  selected  so  that  one  of 
the  mesh  points  determined  from  (20)  did  coincide  with  initially, 
this  would  not  continue  to  hold  as  the  solution  evolved  in  time  causing 
{  }  to  change. 

2)  If  {  }  passes  through  zero,  from  (20)  can  become  arbitrarily  large 

and  the  denominator  of  (20)  becomes  nonanalytic.  Both  of  these  effects 
are  troublesome  numerically. 

3)  Many  of  the  terms  in  I  )  of  (19)  cannot  be  estimated  by  difference 
formulas  involving  values  at  only  two  consecutive  mesh  points.  Since  a 
formula  similar  to  (20)  is  embedded  in  the  differential -algebraic 
system  (13),  the  number  of  consecutive  mesh  points  in  excess  of  two 

used  in  the  formula  increases  the  bandwidth  of  the  Jacobian  and,  consequently, 
the  computation  time. 
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To  circumvent  difficulty  (1),  the  numerator  of  (20)  is  replaced  by  an 
auxiliary  variable  A  which  is  set  initially  and  automatically  adjusted  as  the 
solution  evolves  in  sue h  a  manner  that  the  preselected  number  of  mesh  points 
exactly  fits  into  the  specified  spatial  interval.  With  this  modification  we 


Note  the  change  in  philosophy,  i.e.,  the  number  of  mesh  points  is  selected,  and 
the  spatial  truncation  error  is  a  consequence  of  this  number.  TE  also  varies 
in  time. 

Difficulty  (2)  is  surmounted  by  replacing  the  denominator  of  (20)  by 
F"TT7  where  c  is  a  positive  constant  (typically  0(1)  in  our  numerical 
examples).  This  change  has  the  effect  of  making  the  denominator  analytic  and 
preventing  the  occurrence  of  arbitrarily  large  mesh  intervals  even  if  {  }  -+  0. 

Of  course,  the  smaller  c  is  chosen  the  closer  to  a  nonanalytic  function  v' 
becomes,  and  the  more  numerical  difficulties  can  be  expected. 

For  difficulty  (3)  the  solution  is  less  satisfactory.  Instead  of 
7{  }2  V  c  which  we  cannot  hope  to  evaluate  within  the  constraints  of  our 

formulation  (i.e.,  minimum  bandwidth)  we  employ  a  function  FCx^-  ,v_j  ) 
(already  introduced  in  equation  (8))  which  we  hope  can  be  selected  to  mimic  the 
important  features  of  A  }2+  e.  Selection  of  this  function  is  the  principal 
obstacle  in  the  development  of  an  automatic  adaptive  mesh  method. 


Implicit  in  our  definition  of  the  function  F  is  that  it  does  not  depend  on 
Ti  and  Ti+.|.  This  is  not  an  essential  restriction.  However,  for  the  example 
problems  considered  in  this  study,  T  is  always  in  the  range  [0,1]  and  so  it  was 
felt  that  T  was  not  a  crucial  quantity  in  F.  Further,  spatial  derivatives  of 
T  can  always  be  expressed  in  terms  of  v. 

Another  important  assumption  regarding  F  is  that  it  not  be  a  function  of  any 
of  the  time  derivatives.  The  reason  for  this  is  that  GEARIB  requires  that  all 
time  derivatives  appear  linearly  so  that  the  differential -algebraic  system  can  be 
cast  in  the  form  (13).  The  types  of  F  considered  in  this  study  definitely  do  not 
have  this  property. 
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Of  the  remaining  terms  in  (19)  the  only  one  that  can  be  approximated  in 
terms  of  values  at  only  two  consecutive  mesh  points  is  yvmTn'",  i.e., 


vmT~ 
m  m 


vi  +  Vv 


'vi+i  ~  vi 

v  Hi 


'i+l 


2 

vi 


2Hi 


(22) 


One  could  hope  that  (22)  is  a  reasonable  measure  of  the  discretization  of  the 

2 

nonlinear  term  of  Burger's  equation  even  though  the  other  0(H  )  contribution, 

Tv",  is  neglected, 
mm  3 

In  typical  applications  of  Burger's  equation,  the  coefficient,  S,  of  the 
diffusion  term  on  the  RHS  of  (1)  is  small.  For  our  examples,  it  never  exceeds 
0.01.  Despite  this  small  coefficient,  the  diffusion  term  is  important  at  the 
"knees"  of  the  wave-like  solutions  to  Burger's  equation.  For  example,  see 
Figure  I  of  the  next  section  and  note  the  sharpness  of  the  knees  for  larger 
values  of  t.  As  B  is  reduced  these  knees  become  sharper.  Although  we  are  not 
able  to  estimate  v"'  of  (19)  required  for  the  truncation  error  of  this  diffusion 
term,  we  can  at  least  determine  when  the  diffusion  term  is  important  relative 
to  the  nonlinear  term  by  evaluating  the  ratio  of  these  two  terms,  i.e.. 


p 

Strictly  speaking  the  denominator  of  (23)  should  have  a  factor  (T^  +  T^)  /4, 
but  this  has  been  omitted  for  reasons  cited  previously.  When  R  «  1  the 
diffusion  term  is  not  important  and  (22)  may  be  a  good  measure  of  the  truncation 
error  corresponding  to  the  RHS  of  Burger's  equation.  However,  if  R  _>  0(1), 
then  some  function  should  be  devised  which  forces  an  increased  concentration  of 
mesh  points  within  this  region  where  the  diffusion  term  is  important. 

The  preceding  are  some  of  the  considerations  that  guided  our  selection  of 
the  function  F( x^  ,v .  ,v.+.j )  for  the  determination  of  mesh  spacing.  Numerical 
experiments  were  performed  with  a  variety  of  such  functions.  A  number  of 
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unexpected  phenomena  were  encountered  during  these  experiments,  which  are  described 
in  the  following  section. 

E.  Selection  of  F(x^  >x1+^  ,v.  .v^ )  and  Numerical  Results 

Most  of  the  numerical  results  presented  in  this  section  are  based  on  the 
following  values: 

Y  =  1 ,  3  =  0.01 ,  EPS  =  0.0001 ,  N  =  31 

Y  and  3  are  constant  parameters  in  Burger's  Equation  (1),  EPS  specifies  the  accuracy 
of  the  time  integration  in  GEARIB,  and  N  is  the  number  of  mesh  points.  This  value 
of  EPS  is  small  enough  so  that  the  spatial  discretization  error  dominates  that  of 
the  time  integration. 

In  using  the  GEARIB  integrator,  the  following  were  adhered  to  for  all  of  our 
numerical  work: 

1)  The  Jacobian  of  the  RHS  of  (13)  was  evaluated  analytically  in  the  user 
written  subroutine  PDG.  (There  is  an  option  in  GEARIB  to  evaluate  the 
Jacobian  by  numerical  differencing,  but  the  analytic  approach  is 
generally  preferable.) 

2)  In  computing  the  Newton  matrix,  the  variability  of  the  matrix  on  the  LHS 
of  (13)  was  explicitly  accounted  for.  Some  early  runs  which  did  not  do 
this  suffered  degraded  integration  performance. 

3)  To  start  the  integration,  the  time  derivatives  of  all  dependent  variables 
were  computed  (in  order  to  get  a  good  initial  prediction  and,  consequently, 
improved  corrector  convergence  properties  for  the  first  step.)  Computation 
of  these  initial  derivatives  requires  differentiation  of  all  algebraic 
equations  of  system  (13).  However,  the  coding  for  these  differentiations 
was  already  present  in  subroutine  PDG,  so  relatively  little  additional 
work  was  required. 

Most  of  our  numerical  examples  are  based  on  the  x-interval  [0,2]  with  the 
initial  condition 
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1(1  +  COS (?rX  )  )/2  ,  0  £  X  £  1 

(24) 

0  ,  1  <  x  <  2 

T(x)  has  a  discontinuity  in  its  second  derivative  at  x=l .  Su:h  discontinuities 
sometimes  cause  difficulty  for  nonlinear  iteration  schemes,  especially  in  the 
determination  of  initial  conditions.  However,  derivative  discontinuities  in 
initial  conditions  are  not  uncommon  in  practice,  so  it  was  considered  relevant 
to  use  an  initial  condition  such  as  (24). 

For  the  reasons  discussed  in  Section  C,  the  boundary  conditions  used  are 
always : 


v-|  =0,  =  constant 


(25) 


E.l  Case  I 

Although  not  motivated  by  the  truncation  error  analysis  of  the  preceding 
section,  the  first  numerical  tests  used: 

F^xi  ,xi+l  ,vi  ,vi+l )  ;  \/(vi  +  Vi+1)2/Z  +  E  (26a) 

and 

F(x.,x.+1  ,vi  ,vj+1)  \Avi  +  vi+i )  (26b) 

In  regions  where  v^  is  small  relative  to  ve,  the  mesh  spacing  is  nearly  uniform. 
In  regions  of  steep  slopes  of  the  T-profile,  the  mesh  spacing  is  nearly  inversely 
proportional  to  the  slope. 

Included  in  this  section  are  plots  of  the  numerical  results  for  various 
interesting  and  informative  cases.  We  denote  these  as  Case  I,  Case  II,  ... 

In  Figures  la,  lb,  and  Ic  are  plotted  the  results  for: 

Case  I  N  =  32 

xN  =  1.985 

t  =  1.0 


(27) 
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The  reason  for  the  odd  value  of  x^  is  discussed  following  the  presentation  of 
numerical  results. 

Figure  la  shows  the  T-profiles  at  various  time  intervals  up  to  t  =  2.7, 
which  is  an  arbitrary  stopping  point.  At  this  time  the  influence  of  the  right 
boundary  is  already  strong.  (Later  examples  illustrate  the  effect  of  allowing 
the  integration  to  proceed  to  "indefinitely  large"  values  of  t.)  As  can  be 
seen  from  the  figure,  from  t  =  1.0  to  t  =  2.0  (at  least)  the  wave  moves  to  the 
right  practically  und  storted  and  with  a  constant  speed.  X's  and  +'s  are  used 
on  alternate  profiles  to  indicate  the  solution  points  (x^ ,1^).  In  each  mesh 
interval  an  additional  value  of  T  was  evaluated  at  the  midpoint  using  Hermite 
(cubic)  interpolation  based  on  the  known  quantities  (T^  T^-j,  vi ,  vi+-|).  Each 
T-profile  was  generated  by  connecting  the  N  points  (x^,Ti )  and  the  (N-l) 
intermediate  points  by  straight  lines.  No  symbol  is  used  to  indicate  these 
intermediate  points  on  the  plots.  However,  in  "rough"  regions  of  the  profiles 
these  intermediate  points  are  conspicuous  by  the  presence  of  noticeable  slope 
discontinuities.  These  slope  discontinuities  often  signal  an  inappropriate 
mesh  spacing.  For  example,  in  Case  I  the  mesh  is  too  coarse  in  the  regions  of 
high  curvature  near  the  "knees"  of  the  T  profiles,  and  this  is  clearly  shown  in 
Figure  la.  It  is  apparent  from  the  figure  that  the  x-interval  spanned  by  a 
particular  T-profile  is  less  than  the  total  solution  interval.  The  plot  begins 
at  the  first  x^  for  which 


<IVil  +  'V21’  >  °-2 

and  it  ends  at  the  last  x^  for  which 

( 1  vi- _2 1  +  |vi_1 1)  >  0.2 

This  cutoff  criterion  (determined  experimentally)  is  intended  to  prevent  excessive 
overlapping  of  the  T  profile  plots,  which  makes  it  difficult  to  distinguish 
certain  details. 

Another  type  of  plot  that  was  generated  for  each  case  is  shown  in  Figure  lb. 
It  shows  the  trajectories  of  the  x^  in  time.  This  p’ot  clearly  shows  how  the 
region  of  concentrated  mesh  points  moves  along  with  the  wave.  It  should  be  noted 
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that  a  given  mesh  point  does  not  remain  stationary  with  respect  to  the  wave.  For 
example,  the  trajectory  of  x^g  is  shown  in  Figure  lb  as  a  heavy  line.  Initially 
x^g  lies  in  the  undisturbed  region  ahead  of  the  wave.  As  the  wave  moves  to  the 
right,  x^g  begins  to  move  to  the  right,  but  at  a  speed  less  than  the  wave  speed. 
Consequently,  x^g  passes  smoothly  through  the  wave  region  and  then  slows  nearly  to 
a  constant  asymptotic  position  behind  the  crest. 

It  is  clear  from  this  plot  that  at  any  time  more  than  half  the  mesh  points 
lie  in  a  region  of  nearly  constant  T.  Since  such  regions  of  constant  T  are  not 
of  special  interest  and  do  not  require  many  mesh  points  for  their  description, 
it  is  natural  to  attempt  to  place  more  of  these  mesh  points  within  the  wave  region. 
An  obvious  variant  to  try  is  to  reduce  e  in  (27)  since  the  mesh  size  in  constant 
regions  cannot  exceed  A/«/c.  However,  it  was  found  by  experiment  that  reduction  of 
e  much  below  unity  generally  resulted  in  little  payoff.  The  problem  became  more 
sensitive  numerically  resulting  in  smaller  integration  step  sizes,  more  iterations 
per  step,  and  more  erratic  profiles  and  trajectories.  For  example,  the  slight 
waviness  in  the  x^  trajectories  of  Figure  lb  was  often  accentuated  when  e  was 
reduced. 

A  third  type  of  plot  generated  for  each  case  is  shown  in  Figure  Ic.  It  shows 
the  trajectories  of  the  v^  in  time.  Since  v^  should  be  negative  for  the  examples 
being  considered  in  this  report,  we  have  plotted  -v 1 vimax | ,  where  vimax  is  the 
largest  value  of  v^  (over  all  i  and  all  t)  encountered  during  the  integration. 

From  t  =  1.0  to  t  =  2.0  (approximately)  the  pattern  of  these  v-trajectories  is 
nearly  uniform.  However,  near  t  =  2.7  the  pattern  starts  to  change  because  of  the 
influence  of  the  right  hand  boundary.  Notice  how  v.  corresponding  to  a  position 
ahead  of  the  wave  begins  to  oscillate  about  zero  as  the  wave  approaches,  then  builds 
up  to  its  peak  value  as  its  mesh  point  x.  passes  through  the  wave,  then 
decreases  and  oscillates  about  zero  before  reaching  a  steady  zero  value. 

The  tic  marks  just  above  the  horizontal  axis  of  Figure  Ic  indicate  the 
integration  step  sizes  that  were  used.  To  reduce  the  memory  and  CPU  requirements 
of  the  plotting  phase,  the  trajectories  of  Figures  lb  and  Ic  were  generated  by 
connecting  with  straight  lines  the  values  at  integration  steps  0,  3,  6,  ...  This 
is  especially  conspicuous  in  Figure  Ic.  Several  of  these  straight  line  segments 
spanning  three  integration  steps  are  indicated  by  vertical  lines  at  the  bottom  of 
Figure  Ic. 
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From  the  plot  of  step  sizes,  it  is  apparent  that  starting  at  about  t  =  0.8, 
the  step  size  variation  has  settled  into  a  somewhat  regular  pattern  of  increases 
and  decreases.  There  is  no  obvious  characteristic  of  the  solution  which  would 
require  such  step  size  changes.  More  likely  it  is  a  deficiency  of  the  GEARIB 
step  size  changing  criterion.  This  criterion  could  probably  be  tailored  to  the 
specific  examples  being  considered  in  this  study  in  order  to  reduce  this  wasteful 
variation  of  step  size.  From  t  =  0  to  t  =  0.3,  the  integration  method  was  first 
order.  At  t  =  0.3  GEARIB  switched  (automatically)  to  a  second  order  method  and 
continued  with  that  order  until  the  end.  A  few  runs  similar  to  Case  I  were 
rerun  with  GEARIB  modified  so  that  it  would  always  operate  with  a  first  order 
method.  The  effect  was  that  the  step  size  variation  was  somewhat  less  erratic, 
but  the  average  step  size  was  less,  making  the  integration  less  efficient. 

Beyond  this  experiment  no  further  effort  was  expended  in  tuning  GEARIB  for  the 
class  of  problems  being  studied. 

The  implicit  methods  of  GEARIB  are  implemented  as  an  explicit  prediction 

followed  by  one  or  more  nonlinear  corrections.  Each  of  these  corrections  requires 
the  evaluation  of  the  RHS  of  (13)  and  the  solution  of  a  linear  system.  The 
matrix  of  this  linear  system  is  a  linear  combination  of  the  matrix  on  the  left 
of  (13)  and  the  Jacobian  of  the  RHS.  In  GEARIB  this  matrix  is  re-evaluated  and 
refactored  only  when  the  corrector  convergence  begins  to  deteriorate  badly. 

Thus,  the  number  of  steps  taken  is  not  a  reliable  measure  of  the  efficiency  of 
the  integration.  The  ratios  (NRE/NSTEP)  and  (NJE/NSTEP)  are  also  crucial. 

Here  NRE  is  the  number  of  RHS  evaluations  and  NJE  is  the  number  of  Jacobian 
evaluations.  When  GEARIB  is  operating  efficiently,  typical  values  of  these 
ratios  might  be: 


(fijffp)  =  1.5  to  2.0,  (j^fp)  =  0.05  to  0.1 

In  Case  I  these  ratios  were  4.0  and  0.5,  respectively.  In  some  of  the  other  cases 
to  be  described  later,  these  ratios  improved  but  never  reached  the  hoped-for 
levels. 

Case  I  was  repeated  using  (26b)  and  modifying  (13)  to  use  sums  of  products 
for  terms  such  as  vx  and  vT  (instead  of  products  of  sums).  Th’s  resulted  in  an 
improvement  in  the  smoothness  of  the  T-profiles  and  the  v- trajectories  and  the 
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integration  statistics  were  better.  Based  on  these  results,  this  modification 
of  the  Case  I  formulation  seems  preferable.  However,  the  Case  I  and  "modified 
Case  I"  results  are  based  on  a  somewhat  artificial  example.  In  both  instances, 
was  adjusted  so  that  when  the  initial  x.  were  solved  for  (to  satisfy  (27)), 
one  of  the  interior  x^  fell  exactly  at  x=l .  Then  when  the  initial  v_j  were 
evaluated  based  on  the  obtained  from  (24),  v^  =  0  for  x^  >  1 .  These  smooth 
initial  conditions  were  apparently  the  reason  for  the  superior  performance  of  the 
modified  Case  I.  Later  "more  realistic"  cases  were  run  where  the  right  boundary 
point  was  set  to  2  a  priori.  Subsequent  computation  of  the  initial  x^  resulted 
in  x  =  1  falling  in  the  interior  of  a  mesh  interval  with  the  consequence  that 
v.j+^  =  -  v_j  for  x-  >1.  These  rougher  initial  conditions  had  little  effect 
on  Case  I  but  the  performance  of  the  modified  formulation  was  significantly 
degraded.  Based  on  these  experiments  it  was  felt  that  the  product  of  sums 
formulation  was  preferable.  All  of  the  numerical  results  presented  in  this 
section  are  for  "product  of  sums"  formulations. 


One  interesting  difference  between  Case  I  and  the  modified  formulation  was 
in  the  speed  of  the  wave  as  it  traveled  practically  without  distortion  (say 
between  t  =  1.0  and  t  =  2.0).  Theory  (see  [7])  says  that  the  wave  speed  should 
be  y/2  (=  0.5  in  our  case).  For  Case  I,  the  T-profiles  obey  this  to  as  close  as 
can  be  measured  on  the  plot.  For  the  modified  formulation,  the  wave  speed  is 
0.48.  When  translated  appropriately,  the  profiles  themselves  are  practically 
indistinguishable  in  the  regions  of  steep  slope.  The  reason  for  this  difference 
in  wave  speeds  might  be  expected  to  lie  in  the  discretization  of  the  nonlinear 
term  yvT  on  the  RHS  of  (3a).  For  Case  I  this  is: 


whereas  for  the  modified  Case  I  it  is: 


v<Vi  +  ViW'2  <28b> 

If  (28a)  were  consistently  larger  than  (28b),  this  might  explain  the  larger  wave 
speed  corresponding  to  (28a).  However,  it  is  easy  to  show  that  (28b)  is  larger 
as  long  as  (v^  -  v.+^)  >  0  which  is  over  the  upper  half  of  the  wave,  approximately. 
Thus,  the  reason  for  the  difference  in  wave  speeds  is  not  obvious. 
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E.2  Case  II 

The  oscillations  at  the  crest  and  trough  of  the  T-profiles  of  Figure  1.1 
clearly  indicate  the  inadequacy  of  the  mesh  spacing  formula  (27).  More  mesh 
points  are  needed  in  the  regions  of  high  curvature  at  the  crest  and  trough. 
Probably  fewer  points  are  needed  in  the  steep  portion  of  the  wave.  To  get  this 
increased  concentration  in  regions  of  high  curvature  formula  (22)  was  used  in 
conjunction  with  (19)  (neglecting  all  terms  except  v  T"). 

Case  II 

N  =  31 

*»  '  2-° 
c  =  1.0 


The  initial  T^  again  satisfied  the  formula  (24).  The  boundary  conditions  were 
in  accordance  with  (25).  Figure  I  la  shows  the  T~profiles  for  t  =  0,  0.5,  0.6479. 
The  integration  "hung  up"  at  t  -  0.6479  for  reasons  to  be  discussed  shortly. 
Because  of  a  plotter  malfunction,  the  profiles  between  x  =  0  and  x  =  0.2  are 
distorted.  The  initial  profile  is  identical  to  that  of  Figure  la.  Figures  I I b 
and  lie  show  the  x^  and  v^  trajectories  and  the  integration  step  sizes. 

The  difficulty  with  using  formula  (29)  to  specify  the  mesh  spacing  can  be 
seen  from  Figure  lib.  The  initial  mesh  spacing  seems  to  be  preferable  to  that 
of  Figure  lb,  since  the  mesh  is  more  concentrated  in  regions  where  the  curvature 
of  the  T-profile  is  large.  As  the  wave  moves  to  the  right,  the  mesh  intervals 
(xi+1  ‘  x-j)  corresponding  to  small  i  reach  their  maximum  value  (=  A/‘4»/c~)  and  so 
these  mesh  points  can  no  longer  follow  the  wave.  In  order  to  obtain  the  required 
concentration  of  mesh  points  in  the  region  of  high  curvature  near  the  crest, 
the  following  pattern  occurs: 
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1)  A  mesh  point  moves  out  of  the  region  of  coarse  mesh  in  the  middle  of 
wave  toward  the  crest.  Since  the  lower  part  of  the  wave  is  nearly  a 
mirror  image  of  upper  part,  another  mesh  point  moves  out  of  the  middle 
of  the  wave  toward  the  trough. 

2)  This  continues  until  the  mesh  interval  at  the  center  of  the  wave  reaches 
its  maximum  value  (=  A/^ve").  This  situation  is  characterized  by 

v.  =  V-j+-|,  where  and  x^-j  are  the  mesh  points  in  the  steepest  part 
of  the  wave.  This  is  readily  seen  from  Figure  lie  where  the  v^ 
trajectories  for  the  two  largest  values  of  |v.|  intersect  at  various 
values  of  t.  {The  trajectories  corresponding  to  the  other  large  values 
of  |v.j|  also  intersect  in  pairs  at  the  same  t-points.) 

3)  Now  the  only  way  to  achieve  the  required  mesh  concentration  at  the 
crest  is  for  a  mesh  point  to  cross  over  from  the  trough  to  the  crest. 

As  the  wave  steepens,  these  crossovers  become  more  and  more  violent 
creating  disturbances  throughout  the  field.  These  crossovers  cause 
difficulty  for  the  GEARIB  integrator  as  can  be  seen  from  the  step  sizes 
shown  on  Figure  lie. 

Figure  lib  shows  that  this  pattern  has  been  repeated  four  times.  At  the 
end  of  the  run  (t  =  0.6479)  the  mesh  interval  at  the  center  of  the  wave  has 
reached  its  maximum  and  it  is  time  for  another  crossover  to  begin.  However,  the 
beginning  of  this  crossover  is  accompanied  by  such  violent  changes  in  the 
trajectories  that  the  integrator  continually  failed  to  take  another  successful 
step  even  though  the  step  was  retried  with  ever  smaller  step  sizes.  The  final 
attempt  was  with  a  step  10  orders  of  magnitude  smaller  than  that  of  the  first 
attempt.  Although  it  cannot  be  discerned  from  Figure  lie,  the  GEARIB  integrator 
was  just  barely  able  to  negotiate  the  crossover  beginning  near  t  =  0.42.  The 
final  T-profile  of  Figure  1 1 a  also  illustrates  the  "extreme"  distribution  of 
mesh  points  and  the  v-oscillations  associated  with  this  distribution. 

Several  variations  of  the  mesh  spacing  formula  (29)  were  attempted.  One 

was: 
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This  formula  was  used  with  both  e  =  1.0  and  c  =  10.  After  the  initial  conditions 
converged,  the  integrator  was  unable  to  take  even  one  successful  step.  The  reason 
for  this  appears  to  be  the  following: 

The  initial  conditions  on  v  in  the  region  of  constant  T  (x  =  [1,2])  are  such 
that  v-  =  v^+-|.  The  magnitude  of  these  "spurious"  v..  are  typically  small,  e.g., 

| v . |  =  0(0.05).  However,  the  difference  quotient  in  (30)  based  on  these  vi  is 
0(1),  and  so  A  (from  formula  10)  is  significantly  affected  by  these  small  v- 
oscillations.  Now  the  effect  of  the  diffusive  term  in  Burger's  equation  is  to 
smooth  out  these  high  frequency  oscillations.  Because  there  are  so  many  of  these 
oscillations,  the  cumulative  effect  of  their  being  smoothed  out  is  to  cause  a 
rapid  reduction  in  A,  which  in  turn  forces  a  rapid  change  in  all  the  mesh  points 
through  formula  (30).  The  integrator  is  unable  to  cope  with  these  violent  changes. 

There  are  several  ways  out  of  this  v-oscil lation  difficulty: 

1)  "Doctor"  the  problem  so  that  the  initial  conditions  do  not  contain  these 
spurious  v-oscillations.  For  example,  this  was  accomplished  in  Case  I 

by  adjusting  the  right  boundary,  xN>  so  that  a  mesh  point  fell  exactly  at 
x  =  1.  Such  doctoring  is  not  in  the  spirit  of  developing  a  general 
purpose  numerical  technique  so  nothing  further  was  done  in  this  direction. 

2)  Reduce  the  interval  of  interest  to  x  =  [0,1].  The  initial  conditions 
now  have  no  region  of  constant  T  and  consequently  no  region  of  v- 
oscillation.  This  problem  was  run  with  e  =  10.  There  was  no  longer  any 
difficulty  getting  the  integration  started.  However,  the  integration 
"hung  up"  at  t  =  0.5779  because  of  a  crossover  effect  similar  to  that 
discussed  in  conjunction  with  Case  II.  Due  to  the  shortened  x-interval, 
the  right  boundary  affected  the  wave  from  the  beginning  and  so  the  upper 
and  lower  parts  of  the  wave  were  not  mirror  images  for  t>0.  Consequently, 
the  T-profiles  and  the  x^  and  v^  trajectories  are  somewhat  different  than 
those  of  Figure  II.  The  basic  difficulty  is  still  the  same,  i.e.,  as  the 
lower  part  of  the  wave  straightens  out  (e.g.,  like  the  last  profile  of 
Figure  la)  fewer  mesh  points  are  needed  there,  and  so  from  time  to  time 

a  mesh  point  must  cross  the  steep  portion  of  the  wave  and  into  the  crest 
to  compensate  for  the  mesh  points  that  have  not  been  able  to  follow  the 
wave.  As  the  wave  becomes  steeper,  this  crossover  becomes  more  and  more 
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violent  and  finally  the  integrator  quits.  Of  course,  changing  the 
x-interval  from  [0,2]  to  [0,1]  in  order  to  get  the  integration  started 
is  an  artifice  which  is  in  general  unacceptable. 

3)  A  more  satisfactory  approach  is  to  modify  the  mesh  spacing  formula  to 
nullify  the  effect  of  the  spurious  v-oscillations.  In  fact,  this  has 
already  been  done  in  Case  II  where  the  mesh  spacing  formula  (29)  contains 
the  factor  (v^  +  v^+-j).  Other  variations,  all  containing  this  crucial 
factor,  are  discussed  in  the  following: 

One  mesh  spacing  formula  for  which  a  number  of  tests  were  run  is: 


In  regions  where  T  varies  significantly,  the  ratio  [  ]  in  (31)  is  close  to  unity 

and  so  (31)  is  nearly  equivalent  to  (30).  In  regions  where  T  is  constant, 
vi  =  -vi+1,  and  (31)  is  equivalent  to  (29). 

When  T  is  not  constant  formula  (31)  (in  comparison  to  formula  (29)  of  Case  II) 
has  the  effect  of  putting  a  greater  concentration  of  mesh  points  in  regions  where 
v  is  small  (e.g.,  near  the  crest  and  the  trough)  and  a  lesser  concentration  where 
v  is  large  (e.g.,  in  the  middle  of  the  wave).  The  numerical  results  using  (31) 
were  similar  to  Case  II.  However,  (31)  produced  smoother  T-profiles  and  x  and  v 
trajectories,  and  the  integration  was  about  25%  more  efficient  (for  the  same  value 
of  t).  Equation  (31)  allowed  one  more  crossover  before  becoming  hung  up  (at  t  =  0.86). 
Another  run  using  (31)  with  e  =  10  (instead  of  r.  =  1)  produced  very  similar  results 
but  was  about  25%  more  efficient.  Still  another  run  was  made  with  (31)  and  c  =  10 
but  with  the  number  of  mesh  points,  N  =  51  (instead  of  31).  Crossovers  were  more 
frequent,  but  less  violent.  However,  the  integration  eventually  hung  up  at  t  =  1.05 
(at  the  beginning  of  the  13th  crossover)  for  the  same  reason  as  before.  The 
integration  was  about  25%  less  efficient  (in  terms  of  the  number  of  corrector 
iterations  and  Jacobian  re-evaluations)  than  the  N  =  31,  e  =  10  case.  The  T-profiles 
were  somewhat  smoother  as  would  be  expected  from  the  increased  number  of  mesh 
points. 


1-29 


Bell  Aerospace 


TEXTRON^ 


E.3  Case  III 

In  this  section  we  present  numerical  results  for  a  case  somewhat  similar  to 
the  case  discussed  at  the  end  of  the  preceding  section. 

Case  III 

N  =  51 

xN  ~  2.0 
£  =  10.0 


<xhi  -  xi> 


-  xi^  [«»f  *  vftl 


+  10"30) 


The  difference  between  Case  III  and  that  discussed  at  the  end  of  the  preceding 
section  is  the  use  of  /  in  (32)  instead  of  y  .  The  numerical  results  are 
presented  in  Figures  Ilia,  Illb,  and  IIIc.  The  use  of  /  results  in  quite 
different  looking  x  and  v  trajectories.  There  is  still  the  crossover  effect 
which  occurs  when  the  mesh  spacing  at  the  center  of  the  wave  reaches  its 
maximum  (=  A/v/e).  Now  for  Case  II  and  the  other  similar  cases  discussed  in  the 
preceding  section,  there  was  a  sudden  change  in  the  x  and  v  trajectories  at 
such  maximum  points  (see  Figures  1 1 b  and  lie).  In  contrast,  for  the  present 
case  this  maximum  mesh  spacing  decreases  gradually  and  the  trajectories  are  smooth. 
Then  suddenly  the  mesh  point  crossover  takes  place  and  disturbances  are  propagated 
throughout  the  field.  As  expected,  the  integration  step  size  is  greatly  reduced 
at  such  points. 

From  Figure  1 1 1 b  it  can  be  seen  that  the  integration  hangs  up  at  a  time 
when  the  trajectories  are  smooth,  in  contrast  to  the  cases  previously  discussed, 
e.g..  Figure  lib.  The  reason  that  the  integration  stopped  in  Case  III  was  that 
two  mesh  points  became  equal  to  six  figures  making  the  problem  numerically  too 
sensitive  for  the  integration  to  continue.  The  effect  of  using  f  (instead  of 

4 

/  )  in  the  mesh  spacing  formula  is  to  force  the  mesh  points  even  closer  together 
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Figure  Ilia  -  T-Profiles  for  Case  III 
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in  regions  of  high  curvature.  However,  coincidence  of  mesh  points  to  six  figures 
is  almost  certainly  a  numerical  problem.  One  way  to  avert  such  difficulties  is 
to  modify  (32)  so  as  to  place  a  lower  bound  on  (x-+^  -  x^},  e.g., 

(xi+1  -  x.)  =  A(-~  +  6)  (33) 

/ 

where  6  is  related  to  the  minimum  permissible  mesh  spacing.  This  approach  was  not 
pursued,  however,  since  the  y  formula  seemed  to  be  preferable  (in  terms  of 
integration  efficiency  and  smoothness  of  the  T-profiles)  to  formula  (32)  even 
before  excessively  small  mesh  sizes  became  a  problem. 

E.4  Case  IV 

Problems  arose  in  the  examples  of  the  preceding  subsections  due  to  the  fact 
that  as  the  wave  moved  to  the  right,  additional  mesh  points  were  needed  behind 
the  wave.  These  mesh  points  had  to  be  supplied  from  those  initially  ahead  of  the 
wave,  and  the  only  way  they  could  get  to  where  they  were  needed  was  to  pass 
through  the  wave.  As  the  wave  became  steeper,  this  passage  became  more  and  more 
difficult  and  finally  prevented  the  integration  from  continuing. 

This  suggests  adding  mesh  points  from  the  left  boundary,  as  required,  making 
it  unnecessary  for  mesh  points  to  pass  through  the  wave.  To  do  this  "properly" 
would  require: 

1)  Establishing  the  criterion  by  which  a  new  mesh  point  could  be  admitted 
from  the  left  boundary.  Most  likely  this  would  be  based  on  how  near 
(x£  -  x-j )  is  to  its  maximum  allowable  value.  It  would  be  necessary  to 
monitor  for  such  a  condition  frequently. 

2)  Since  (x^+-j  -  x^)  appears  in  the  denominator  of  some  of  the  terms  of  the 
governing  di f ferential -al gebraic  system  (13)  it  would  be  necessary  to 
develop  special  formulas  to  allow  the  newly  introduced  mesh  interval  to 
start  off  with  a  zero  value.  Alternatively,  the  new  mesh  interval  might 
start  off  with  a  finite  value  (say  by  halving  the  interval  (x^  -  x^ 
that  existed  just  prior  to  the  admission  of  the  new  mesh  point). 

However,  this  would  require  the  specification  of  new  values  of  v  and  T 

at  this  new  interior  field  point  such  that  all  of  the  algebraic  equations 
are  satisfied. 
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3)  Shift  all  the  xi ,  ,  Ti  (i  =  2,3,...,N)  and  A  of  the  dependent  variable 

vector,  add  the  x,  v,  and  T  values  for  the  new  mesh  point,  and  then 
restart  the  integration. 

Although  this  approach  is  probably  possible,  it  would  require  a  considerable  amount 
of  analysis  and  reprogramming  to  accomplish  it.  A  simple  alternative  approach  is 
to  force  x-|  to  move  to  the  right  in  time.  This  is  easily  accomplished  by  prescribing 
f, (t)  on  the  right  of  the  second  equation  of  (13)  to  be  a  positive  monotonic 
function  of  t.  As  noted  previously,  theory  says  that  the  wave-like  solutions  of 
Burger's  equation  travel  at  a  speed  y/2  (=  0.5  in  the  present  case).  We  then 
define: 

Case  IV 

N  =  31 

x-j  =  t/2 

XN  =  2,0 
c  =  1.0 


Note  that  (34)  is  identical  to  formula  (31)  discussed  in  Section  E.2.  The 
numerical  results  are  illustrated  in  Figures  IVa,  IVb,  and  IVc.  From  Figure  IVb, 
it  can  be  seen  that  "crossovers"  still  occur.  However,  the  presence  of  the  moving 
left  boundary  reduces  their  frequency.  When  the  identical  case  was  run  with  a 
fixed  left  boundary,  the  integration  "hung  up"  at  t  =  0.86  (near  the  end  of  the 
fifth  crossover).  The  present  case  hung  up  at  t  =  1.10  (near  the  end  of  the  third 
crossover).  Further,  the  integration  was  about  50%  more  efficient  (for  the  same 
values  of  t)  for  the  present  case. 


1-34 


Bell  Aerospace 


TEXTRON 


E.5  Case  V 

Forcing  the  left  boundary  to  move  with  the  wave  speed  reduced  the  crossover 
effect  and  allowed  the  integration  to  proceed  further.  However,  the  problem  is 
still  essentially  the  same  as  before,  i.e.,  as  the  wave  moves  to  the  right  there 
is  an  excessive  number  of  mesh  points  ahead  of  the  wave  (relative  to  those  behind 
the  wave).  An  obvious  remedy  for  this  is  to  remove  mesh  points  from  the  right 
boundary  as  they  are  no  longer  needed.  Such  removal  of  mesh  points  entails  some 
of  the  same  difficulties  cited  in  Section  E.4  when  mesh  points  are  added  through 
the  left  boundary.  Generally  though,  removal  of  mesh  points  is  somewhat  easier 
than  addition.  This  is  discussed  further  in  the  next  section.  For  the  present, 
we  simulate  the  removal  of  mesh  points  by  forcing  the  right  boundary  to  move 
with  the  wave  speed. 

Case  V 

N  =  31 
x1  =  t/2 
xN  =  2.0  +  t/2 

E  =  1.0 

(xi+l  -  x.):  Formula  (34),  Same  as  Class  IV 

The  numerical  results  are  shown  in  Figures  Va,  Vb,  and  Vc.  Because  of  a  plotter 
malfunction,  the  beginning  (near  x  =  0)  of  Figure  Va  is  slightly  distorted.  In 
Figures  Vb  and  Vc,  the  plots  end  near  t  =  2.7  even  though  the  integration  went 
to  at  least  t  =  2.9.  This  is  due  to  a  plot  program  idiosyncracy  which  does  not 
plot  the  last  integration  point,  which  was  reached  by  a  large  integration  step 
(At  >  0.2).  Also,  the  uniform  region  above  x  =  2.2  on  Figure  Vb  has  been  trimmed 
off. 

From  the  plots  it  is  apparent  that  the  use  of  boundary  points  which  move 
at  the  wave  speed  has  eliminated  problems  associated  with  crossovers.  The  mesh 
now  moves  along  with  the  wave  and  all  dependent  variables  behave  either  as 
constants  or  linearly  in  time  allowing  large  integration  steps  to  be  taken.  The 
large  reduction  in  step  size  near  t  -  1.10  does  not  appear  to  be  related  to  any 
properties  of  the  solution  and  so  should  be  attributed  to  a  numerical  integration  fluki 
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E.6  Case  VI 

Although  the  use  of  moving  boundaries  in  Case  V  allowed  the  integration  to 
proceed  smoothly,  it  is  not  a  satisfactory  solution  for  the  following  reasons: 

1)  Knowledge  of  the  direction  and  speed  of  the  wave  is  required.  In  more 
general  problems,  e.g.,  gas  dynamics,  this  information  would  not 
generally  be  known. 

2)  The  boundary  conditions  on  v  and  T  are  applied  at  the  moving  boundaries 
instead  of  at  the  original  fixed  boundaries.  Thus,  the  problem  actually 
solved  is  not  the  problem  originally  posed.  For  Case  V,  there  is 
probably  little  difference  between  the  actual  problem  and  the  problem 
solved  until  the  wave  begins  to  be  affected  by  the  downstream  boundary 
(at  about  t  =  2.5  according  to  Figure  la). 

To  overcome  objections  (1)  and  (2)  but  still  retain  the  desirable  properties 
of  Case  V,  the  following  modification  of  system  (13)  was  considered: 

1)  For  the  first  few  mesh  intervals  (say  l,2,...,m)  use  one  of  the  "standard1 
types  of  formulas,  e.g.,  (29),  (30),  (31),  but  with  a  reduced  value  of  e 
(say  e).  To  indicate  this  we  rewrite  the  general  mesh  spacing  formula 
(8)  as: 


(xi+l  -  xi>  ' 


The  examples  so  far  have  used  an  e  =  1 .  If  an  e  =  1/16  is  used  in  an  F 
containing  V  then  the  maximum  possible  mesh  size  is  double  what  it  is 
with  e  =  1.  We  have  previously  mentioned  some  numerical  experiments  for 
which  small  values  of  e  have  been  used  for  all  intervals.  These  experiments 
were  generally  unsuccessful.  However,  it  was  hoped  that  by  using  a  reduced 
e  in  only  a  few  intervals  behind  the  wave  the  difficulties  associated  with 
small  e  would  be  minimal.  This  modification  allows  mesh  points  behind  the 
wave  to  follow  the  wave  more  closely  and  so  is  analogous  to  Case  IV 
(moving  left  boundary).  Recall  (Figure  IVb)  that  this  only  delayed 
integration  "hang  up"  due  to  the  crossover  effect.  Also  needed  is  a 
technique  by  which  mesh  points  can  be  removed  from  the  right  boundary. 
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2)  To  permit  mesh  points  ahead  of  the  wave  to  move  toward  and  out  of  the 
right  boundary,  the  usual  mesh  spacing  formula  for  the  last  interval 
[xN_r  xN]  was  replaced  by: 

(x2  -  Xl)  +  (x3  -  x2)  +  ...  +  (xm+1  -  xm)  +  (xN  -  xN_n )  = 


1 


+  ... 


|  F(X]  ,x2,v-j  ,v2»e)  F(x2,x3,v2,v3,e) 

+  F(xm,Vi,vm,Vis£)  F(xN-rVVrVc)  (  (36) 


Note  that  the  standard  c  (rather  than  the  reduced  e)  is  used  in  (36).  In 
our  examples  we  have  used  e  =  1  and  e  =  1/16.  Consider  the  special  case  where 
m  =  1  in  (36).  As  v-|  and  v2+0,(x2  -  x-j  )-*-2A  (from  (35)),  and  the  contribution 
from  the  first  term  on  the  RHS  of  (36)  approaches  A.  The  last  term  on  the  RHS 
of  (36)  cannot  exceed  A.  Thus,  in  order  that  (36)  be  satisfied  (x^  -  x^ _^)  must 
go  to  zero  (or  even  negative  values  if  the  last  term  on  the  RHS  is  less  than  A). 

A  numerical  test  of  this  case  was  run.  (x^  *  ])  approached  zero  asymptotically. 

Because  (xN  -  x^)  appears  in  the  denominator  of  certain  terms  of  the  differential 
algebraic  system  (13),  the  integration  hung  up  because  of  this  ill  conditioning. 


If  m>l  then  (x^  -  x^_j)  must  eventually  become  negative  as  v^ ,  v2,  ...,  v i 
become  small.  There  may  be  some  question  as  to  the  validity  of  the  solution  of 
system  (13)  (modified  by  (35)  and  (36))  for  negative  values  of  (x^  -  xN  -| ) .  This 
did  not  seem  to  be  a  problem  in  the  examples  attempted  with  this  formulation. 
However,  the  integration  for  all  these  examples  hung  up  before  the  wave  approached 
the  right  boundary.  After  presentation  of  numerical  results  for  a  typical  case  of 
this  modified  formulation,  we  will  discuss  a  correct  method  of  treating  vanishing 
mesh  intervals  at  the  right  boundary. 


♦Because  (36)  directly  couples  the  dependent  variable  at  many  mesh  points  (instead  of 
just  two)  it  is  necessary  to  place  this  relation  in  the  last  row  of  (13).  The  last 
row  of  the  Jacobian  of  (13)  will  consequently  have  many  no  zero  elements.  However, 
as  previously  noted  special  variants  of  the  GEARIB  routines  have  been  developed  to 
accommodate  these  cases. 
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Case  VI 

System  (13)  modified  by  (35)  and  (36) 
N  =  31 
x-|  =  0 

XN  =  2,0 

e  =  i.O 

e  =  1/16 

m  =  5 


(xi+i  -  n> 


(37) 


*  W‘ 


4 


h*  >  -  vff 

(vi  +  vi+/ 

\xi+l  "  xi / 

2(v?  +  v^+1  +  IQ"30) 

+  e 


Of  course  (37)  is  modified  in  accordance  with  (35)  for  mesh  intervals  1  to  5  and 
(36)  is  used  for  the  last  mesh  interval.  The  numerical  results  are  shown  in 
Figures  Via,  b,  and  c.  Note  the  plotter  malfunctions  which  caused  the  distortion 
at  the  beginning  of  Figure  Via  and  the  skip  near  the  end  of  Figures  VIb  and  Vic. 

From  Figure  VIb  it  can  be  seen  that  the  trajectories  for  xN-1  and  xN_2  have 
crossed  x^  (=  2.0)  before  the  integration  hung  up.  Although  there  is  a  potential 
danger  of  the  integration  hanging  up  (due  to  a  singularity  in  the  governing  system) 
when  x^_i  =  x^,  this  did  not  occur  for  this  case  because  the  integration  step 
size  used  at  this  point  was  such  that  |x^  -  x^-j)  did  not  become  excessively 
small.  Comparing  Figure  VIb  with  Figure  Vb,  we  see  that  at  least  initially  (up  to 
about  t  =  0.40)  our  new  formulation  has  the  desired  effect.  Near  t  =  0.40, (X2  -  x-j ) 
is  close  to  its  maximum  permissible  value  (=  2A)  but  then  a  slight  spurious  build 
up  of  v 2  (from  0.014  to  0.056)  causes  a  dip  in  the  x2  trajectory.  (This  small 
change  in  v2  has  a  noticeable  effect  on  (x2  -  x-j)  because  of  the  small  value  of 
e  (=  1/16)  being  used  in  (37).)  The  integration  proceeds  with  both  (x2  -  x^)  and 
(x^  -  x2)  near  their  maximum  permissible  values  at  about  t  =  0.61.  In  order  for 
the  mesh  to  continue  to  follow  the  wave,  it  is  now  necessary  for  (x^  -  x^)  to 
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increase.  However,  has  not  fallen  off  enough  to  allow  this.  The  only  way  to 
obtain  the  required  mesh  concentration  in  the  crest  is  for  a  crossover  to  occur. 

From  Figure  VIb,  it  may  not  seem  that  a  crossover  should  be  occurring  since 

the  mesh  size  in  the  middle  of  the  wave  is  nowhere  near  its  maximum  permissible 

size.  However,  recall  that  the  mesh  spacing  formula  (37)  contains  the  additional 
2 

term  (v^  +  vi+-j)  .  There  is  no  way  this  term  can  become  zero  in  the  middle  of 

the  wave.  What  has  happened  is  that  v.  has  become  equal  to  v.+-j  (see  Figure  Vic) 

reducing  the  second  term  in  of  (37)  to  zero.  This  is  the  same  effect  that  we 

have  already  seen  in  Cases  II,  III,  and  IV.  Although  the  modified  formulation  has 
succeeded  in  delaying  the  first  major  crossover,  when  this  crossover  is  required 
it  is  so  violent  that  the  integration  cannot  continue. 

To  eliminate  the  possibility  that  our  questionable  treatment  of  the  problem 
(i.e.,  allowing  the  mesh  size  to  become  negative  for  the  last  interval)  influenced 
the  preceding  conclusion,  the  same  problem  was  solved  in  a  different  manner. 

Instead  of  allowing  to  cross  x^,  the  integration  continued  until  x^  =  xN* 
(to  within  the  user  specified  integration  accuracy).  Then  T^,  vN,  and  xN  were 
eliminated  from  the  system,  and  the  integration  was  restarted  with  the  boundary 

conditions  now  on  T^_-j  and  x^  -j  and  with  a  formula  similar  to  (37)  now  specifying 

the  mesh  spacing  for  the  interval  (x^  -  x^_2).  The  integration  then  continued 
until  x^_2  =  x^_i  where  another  restart  was  made,  etc.  Until  x^  -j  =  x^,  the 
numerical  results  were  identical  to  those  of  Case  VI.  Upon  restart,  there  was  a 
slight  change  in  slopes  of  the  xi  and  v^  trajectories,  especially  for  the  larger 
values  of  i .  At  t  =  0.577,  x^_2  =  x^_-| ,  whereas  in  Case  VI  this  did  not  occur 
until  t  =  0.631.  The  integration  hung  up  at  t  =  0.659  versus  t  =  0.658  in  Case  VI. 
The  T-profiles  and  the  x^  and  v.  trajectories  (except  for  i  close  to  N)  were 
practically  indistinguishable  for  the  two  cases. 


*Near  such  singular  points  the  integrator  has  considerable  difficulty.  The  technique 
by  which  such  singular  points  are  determined  could  be  improved  by  making  use  of  the 
polynomial  representation  (in  t)  of  each  of  the  dependent  variables  which  is 
generated  by  GtARIB.  However,  this  was  not  our  objective  in  this  example  so  no 
effort  was  made  in  this  direction. 
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Several  other  variations  of  this  formulation  were  attempted.  Instead  of  mesh 
spacing  formula  (37),  formula  (34)  was  used.  This  did  not  make  any  significant 
difference.  Other  cases  were  run  using  l  =  1/256  (instead  of  1/16)  in  conjunction 
with  both  formulas  (34)  and  (37).  The  integration  managed  to  proceed  somewhat 
farther  (e.g.,  t  =  0.80)  but  (xg  -  x-j)  was  more  erratic,  as  expected,  because  of 
the  small  value  of  e. 

E.7  Case  VII 

Of  all  the  mesh  spacing  formulas  discussed  so  far  that  of  Case  I  (equation  (27)) 

produced  the  best  results.  However,  as  discussed  previously,  this  formula  is 

clearly  inadequate  near  the  knees  of  the  T-profile.  It  is  almost  certain  that  we 

2  2 

need  a  term  containing  (v.^  -  )  /(x.+-j  -  x^)  ,  but  so  far  all  of  our  attempts 
to  use  this  term  (except  for  the  "artificial"  Case  V)  have  been  thwarted  by  the 

"crossover  effect."  This  suggests  combining  the  term  (v.  +  v,.,,)  of  formula  (27) 

2  2  ‘  1  + 
with  (v.^,  -  v.)  /lx..-,  -  x.)  to  yield  a  mesh  spacing  formula  with  the  desirable 

properties  of  Case  I  (i.e.,  allowing  mesh  points  to  move  smoothly  through  the  wave) 

with  those  of  Cases  II  -  VI  (i.e.,  concentration  of  mesh  points  at  the  knees  of  the 

wave) . 

This  has  already  been  attempted  in  formula  (37)  of  Case  VI.  One  clear  effect 

p 

of  the  additional  term  (v^-  +  v.^ )  in  (37)  is  to  force  a  more  concentrated  mesh 
in  the  middle  of  the  wave.  (Compare  Figure  VIb  with  Figures  1 1 b  or  IVb.)  However, 
this  increased  mesh  concentration  has  not  allowed  mesh  points  to  pass  smoothly 
through  the  wave.  The  integration  in  Case  VI  was  still  forced  to  stop  because  of 
the  violent  crossover  effect  associated  with  (vi+-|  -  )  /(x^  -  x^ )  .  An 

examination  of  the  relative  contributions  of  the  terms  (37)  showed  that  as  the  wave 

2  2  2 
steepened  the  (v.,,  -  v.)  / ( x  - . 1  -  x.)  term  generally  dominated  the  (v.  +  v.,,) 

i+ii  i+ii  2  ii+i 

term,  even  in  the  middle  of  the  wave  where  (v^  +  v^r  is  at  its  largest.  The 
consequence  is  that  the  "easiest"  way  to  provide  the  required  increased  mesh 
concentration  at  the  crest  and  trough  is  to  maximize  the  mesh  interval  at  the  center 
of  the  wave  by  making  v..  =  v..+-j,  which  is  what  happens  at  the  end  of  the  integration 
in  Case  VI  (see  Figure  Vic)  and  also  in  Cases  II,  III,  and  IV. 
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It  seems,  then,  that  to  prevent  the  crossover  effect  the  contribution  of 
2 

(vi  +  V-j.i)  should  be  such  that  this  term  is  dominant  in  the  middle  of  the  wave. 

i  i  + 1  2  p 

The  formula  should  also  be  such  that  (v.^  -  v.)  /(x^  -  x. )  is  important  at  the 
knees  of  the  wave.  This  suggests  the  following  mesh  spacing  formula: 


Here  y  and  8  are  the  coefficients  in  Burger's  equation  (3a),  c,  is  a  constant  to 

be  set  so  that  [  ]  does  not  come  too  close  to  becoming  singular,  and  C  is  a 

numerical  "tuning"  constant.  If  C  =  1  and  =  0,  the  [  ]  term  approximates 

the  square  of  the  ratio  of  the  diffusive  term  to  the  nonlinear  term  of  Burger's 

2 

equation.  (Strictly  speaking  there  should  be  a  factor  (T^  +  w  /A  in  the 

denominator  of  [  ].)  The  rationale  for  (38)  is  that  if  the  diffusive  term  is 

2  2 

important  (e.g.,  at  the  knees  of  the  wave)  then  (vi+1  -  v^ )  /(xi+^  *  )  should 

have  a  significant  influence  in  the  mesh  spacing  formula.  Accordingly,  [  ] 

has  been  defined  so  that  [  ]  >  0(1)  when  the  diffusive  term  is  important. 


Case  VII 

=  31 

=  0. 

=  2.0 
=  0.5 
=  0.1 
=  0.25 
(xi+l  *  xi ) :  formula  (38) 

The  numerical  results  are  shown  in  Figures  Vila,  b,  and  c.  Because  of  plot  routine 
storage  limitations,  the  plots  of  the  x  and  v  trajectories  were  terminated  a  few 
steps  before  the  final  station,  t  =  3.1.  From  the  T-profiles  it  can  be  seen  that 
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Figure  Vila  -  T-Profiles  for  Case  VII 


Figure  V 1 1 c  -  v-Trajectories  and  Integration  Step  Sizes  for  Case  VII 
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t  =  3.1  is  nearly  steady  state.  A  later  example  (Case  IX)  will  show  the  actual 
attainment  of  steady  state. 

A  comparison  of  Figures  Vila  and  b  with  Figures  la  and  b  shows  that  formula 
(38)  does  in  fact  have  the  desired  effect.  The  "crossover  effect"  has  been 
eliminated,  i.e.,  the  mesh  points  are  now  able  to  pass  smoothly  through  the  wave. 
There  is  also  an  increased  concentration  of  mesh  points  at  the  knees  of  the  wave 
which  eliminates  most  of  the  spurious  oscillations  in  the  T-profiles  that  are 
present  in  Figure  la.  The  v-trajectories  of  Figure  VI Ic  look  much  different  than 
those  of  Figure  Ic.  This  is  due  to  the  fact  the  integration  proceeded  slightly 
further  in  Case  VII  (t  =  3.1  versus  t  =  2.7).  The  T-profiles  steepened 
significantly  between  t  -  2.7  and  t  =  3.1.  Consequently,  the  scaling  factor  for 
the  v-trajectories  was  quite  different  in  the  two  cases. 

A  check  of  the  relative  importance  of  the  terms  in  (38)  shows  that  when  the 

2 

wave  is  traveling  undistorted  (from  t  =  1.0  to  t  =  2.5)  the  (v.  +  v.  , )  /4  term 

2p  1 

dominates  the  [  ]  { -  v^  /(xi+1  -  )  term  in  the  mesh  interval  at  the 

"center"  of  the  wave  by  3:1.  However,  one  or  two  mesh  intervals  to  either  side 
of  this  center  interval  the  dominance  is  reversed. 


Comparing  the  integration  statistics  for  Case  I  and  Case  VII  at  t  =  2.7,  we 
find  that  Case  VII  is  approximately  15%  less  efficient  in  terms  of  the  number  of 
corrector  iterations  required.  Case  VII  is  slightly  more  efficient  in  terms  of 
the  number  of  Jacobian  evaluations  required. 

Several  variations  of  Case  VII  were  run.  In  one  set  of  experiments,  the 
effect  of  e.|  was  examined.  Values  of  1.0,  0.1,  and  0.01  were  used.  The  integration 
statistics  for  =  1.0  were  significantly  worse  than  for  e-j  =  0.1.  In  the  case 
of  =  0.01,  the  integration  hung  up  early  in  the  run.  Based  on  these  results 
e-|  =  0.1  was  used  in  all  subsequent  runs. 

Another  set  of  experiments  tested  the  effect  of  varying  e.  e  =  0.5  was 
found  to  be  significantly  better  than  e  =  1.  An  attempt  to  use  e  =  0.25  failed 
to  converge  on  the  initial  conditions.  In  another  test,  an  e  that  decayed  in 
time  was  used  (e  =  0.5e_t).  Initially,  this  was  identical  to  Case  VII.  At  t  =  0.5 
the  variable  e  case  was  statistically  slightly  better  than  Case  VII.  However,  by  the 
end  of  the  run,  its  statistics  were  much  worse  than  those  of  Case  VII.  The  small 
value  of  e  at  the  end  of  the  run  did  have  the  effect  of  forcing  more  mesh  points 
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into  the  wave.  For  example,  the  T-profile  at  t  =  3.0  of  Figure  Vila  shows  11  mesh 
points  in  the  region  x  1.8,  whereas  the  small  e  case  has  17  in  the  same  region. 

Another  set  of  experiments  examined  the  effect  of  using  smaller  values  of  B. 

Reducing  B  makes  the  T-profiles  steeper  and  generally  makes  the  problem  more 

difficult.  A  repeat  of  Case  VII  with  $  =  0.001  (instead  of  0.01)  integrated  very 

inefficiently  and  the  T-profiles  had  significant  oscillations.  This  case  was  rerun 

with  the  constant  C  in  (38)  increased  by  100  to  compensate  for  the  smaller  value 

2 

of  B.  (The  product  CB  is  now  the  same  as  for  Case  VII.)  Up  to  about  t  =  0.1, 
this  case  behaved  similar  to  Case  VII.  Beyond  this  the  integration  statistics 
degraded.  The  integration  was  stopped  (deliberately)  at  t  =  0.6.  By  this  time, 
the  T-profile  was  significantly  steeper  than  for  B  =  0.01.  (Ma  'mum  slope  about 
three  times  as  great.)  The  T-profile  also  had  noticeable  oscillations  at  t  =  0.6 

In  an  attempt  to  improve  the  numerical  results  for  the  B=  0.001  case,  this 
case  was  rerun  using  a  variant  of  the  mesh  spacing  formula  (38)  in  which  the 
exponent  of  the  (v^  +  v^)  term  in  the  denominator  of  [  ]  was  4  (instead  of  2). 

This  modification  reduced  the  spurious  oscillations  in  the  T-profiles.  However, 
the  integration  statistics  at  t  =  0.6  (the  end  of  the  run)  were  about  the  same 
as  they  were  for  an  exponent  of  2.  This  modified  formula  was  also  used  for  a 
B  =  0.0001  case.  The  integration  hung  up  at  t  =  0.51.  Prior  to  t  =  0.51,  the 
T-profiles  and  x  and  v-trajectories  were  quite  erratic.  A  rerun  of  this  B  =  0.0001 
case  using  (v^  +  v^+-|)^  in  [  ]  of  (38)  produced  even  worse  results. 

E.8  Case  VIII 

In  this  section  we  describe  the  numerical  results  for  an  example  problem  for 
which  an  exact  analytical  solution  is  known.  It  is  taken  from  Sincovec  and  Madsen 
[8],  who  also  solved  this  problem  numerically.  This  exact  solution  can  be  written: 


where 


Tex(t’x>  ’ 


(0.1e'A  »  0.5e'g  *  e~c) 
(e  +  e  +  e  ) 


(39) 


A  = 
B  s 


C 


-  (x  -  0.5  +  4.95t)/203 

-  (x  -  0.5  +  0.75t)/4B 

-  (x  -  0.375J/2B 


0  is  the  coefficient  of  the  diffusion  term  in  Burger's  equation  (1).  The  coefficient, 
Y,  of  the  nonlinear  term  is  assumed  to  be  unity  in  (39). 

For  our  examples,  the  x-interval  is  taken  to  be  [0,1]  and  the  boundary  conditions 

are: 

V1  =  0-  W'-1*  (40) 

The  initial  condition  on  T  is  (39)  evaluated  at  t  =  0.  The  initial  conditions  for 
x  and  v  are  computed  in  terms  of  this  analytical  representation  of  Tex(0»x)  just  as 
they  were  previously  for  the  initial  condition  (24). 

Strictly  speaking  the  boundary  condition  v^  =  0  is  not  consistent  with  (39). 
However,  for  the  values  of  3  being  considered,  [3T  /3x](t,0)  from  (39)  is  less 
than  10  .  Initially,  Tex ( t , 1 )  =  0.1  to  eight  figures.  However,  as  the  wave  moves 

to  the  right  this  is  no  longer  true.  Consequently,  the  exact  time  dependent 
boundary  condition  (40)  was  used. 
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The  numerical  results  are  shown  in  Figures  Villa,  b,  and  c.  Because  of  a  flaw  in 
the  logic  that  determines  the  scaling  factor  for  the  v-trajectory  plots,  these 
trajectories  have  not  been  normalized  so  that  the  maximum  value  is  unity.  Also, 
because  of  memory  limitations  for  the  storage  of  plotter  data,  only  every  fourth 
integration  point  has  been  plotted  in  Figures  VUIb  and  VIIIc. 

Note  that  the  initial  T  profile  in  Figure  Villa  has  regions  of  steep  slope  and 
high  curvature,  in  contrast  to  the  gentler  initial  T-profile  used  in  previous 
examples.  Decreasing  3  accentuates  the  slope  and  curvature  of  this  initial  profile. 
We  had  intended  to  use  3  =  0.003  in  Case  VIII  in  order  to  make  a  direct  comparison 
of  our  numerical  results  with  those  of  [8]  (for  which  6  =  0.003  was  used).  It 
turned  out,  however,  that  it  was  difficult  to  get  the  iterations  for  the  initial 
x.j  and  v.j  to  converge  for  such  small  values  of  3.  By  adjusting  the  constants  e-j 
and  C  in  the  mesh  spacing  definition  and  by  adding  the  "lower  bound"  term,  (0.05A), 
in  (41)  and  by  making  some  other  adjustments  in  the  initial  condition  iteration 
scheme,  convergence  could  usually  be  attained.  However,  for  3  =  0.003,  the  initial 
distribution  of  mesh  points  was  erratic  and  the  initial  vi  contained  some  rather 
large,  0(0.3),  spurious  oscillations.  This  was  true  in  spite  of  the  fact  that  the 
residuals  (of  all  the  algebraic  equations  that  must  be  satisfied  by  the  initial 
values  of  T^ ,  v^ ,  xi  and  A)  were  driven  to  very  small  values.  The  integration  was 
usually  able  to  get  started  in  such  cases  but  always  hung  up  for  small  t  as  the  v^ 
oscillations  built  up,  instead  of  dying  off  as  they  normally  do.  In  the  many 
attempts  that  were  made  to  run  this  exact  solution  example  for  small  values  of  3, 

3  =  0.005  was  the  smallest  value  for  which  reasonable  initial  conditions  and 
integration  efficiency  could  be  obtained. 

The  T-profiles  of  Figure  Villa  are  smooth  with  no  significant  spurious 
oscillations  and  an  apparently  good  distribution  of  mesh  points.  The  final  T-profile 
(at  t  =  1.5)  was  the  constant  steady  state  value,  T=1 ,  and  was  not  plotted.  The 
agreement  between  these  profiles  and  the  exact  solution  was  good.  Some  of  the 
largest  discrepancies  occurred  for  the  profile  at  t=l  : 


^Numerical 

TExact 

i 

i 


0.6072 

0.3663 

0.2574 
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For  these  values  of  T  the  slope  is  large  (0(15))  and  the  mesh  spacing  is  small 
(Ax  =  0(0.01)).  Consequently,  a  plot  of  the  exact  solution  would  hardly  be 
distinguishable  from  the  profiles  of  Figure  Villa. 

The  plots  of  the  x  and  v-trajectories  seem  to  be  well  behaved  except  possibly 
for  the  x-trajectories  for  small  t  and  for  1.0  <  t  <  1.2.  Shortly  after  t  =  1.2 
the  steady  solution  has  been  reached. 

From  the  plot  of  step  sizes  of  Figure  VIIIc,  it  is  clear  that  the  integration 
did  not  proceed  efficiently.  Not  only  were  many  steps  taken  but  the  ratio  of 
corrector  iterations  to  steps  was  high  (NRE/NSTEP=»4)  and  the  ratio  of  Jacobian 
evaluations  to  steps  was  high  (NJE/NSTEP-0.5) .  It  is  likely  that  these  statistics 
could  be  improved  somewhat  by  a  better  selection  of  the  tuning  parameters,  e,  , 
and  C.  In  Case  VIII,  these  parameters  were  chosen  to  provide  good  convergence 
properties  of  the  iteration  for  initial  conditions.  By  smoothly  altering  them 
during  the  course  of  the  integration  (toward  the  values  used  in  Case  VII,  say)  a 
more  efficient  integration  may  have  resulted. 

The  numerical  results  of  reference  [8]  were  generated  by  discretizing  the 

standard  Burger's  equation  (1)  on  a  uniform  mesh  (200  points)  using  three  point 

differencing.  The  time  integration  was  performed  by  GEARB  ,  a  subroutine  package 

similar  to  GEARIB.  The  numerical  examples  of  [8]  were  for  B  =  0.003  and  an 

-5  -4 

integration  error  tolerance  of  10  .  (We  used  10  .)  These  two  differences  make 

the  numerical  problem  of  [8]  somewhat  more  demanding  with  respect  to  computation 
time. 

The  time  integration  in  [8]  was  carried  out  to  t  =  1.1.  A  comparison  of 
integration  statistics  with  Case  VIII  at  this  time  is: 

_ NSTEP _ NRE _ NJE 

[8]  434  648  6 

Case  VIII  532  1909  268 


The  maximum  error  in  T  for  the  two  cases  was  about  the  same  (0.02).  The  minimum 
mesh  size  in  Case  VIII  was  about  0.007;  slightly  larger  than  the  uniform  0.005 
mesh  size  of  [8].  (The  0.05  term  in  the  mesh  spacing  formula  (41)  may  have 
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influenced  this  minimum  mesh  size  even  though  the  minimum  allowable  mesh  size  was 
0.05  A  »  0.0023.) 


Now  [8]  uses  many  more  mesh  points  (200)  than  Case  VIII  which  uses  only  31. 
However,  Case  VIII  has  three  dependent  variables  per  mesh  point,  whereas  [8]  only 
has  one.  Thus,  the  ratio  of  the  system  sizes  is  about  2:1.  A  rough  estimate  then 
is  that  it  takes  [8]  about  twice  as  long  to  perform  a  single  corrector  iteration. 
On  this  basis  alone  [8]  seems  to  have  a  definite  advantage  in  terms  of  integration 
efficiency.  Add  to  this  the  great  disparity  in  the  number  of  Jacobian  evaluations 
and  the  fact  that  the  problem  solved  in  [8]  was  more  demanding  computationally. 

The  conclusion  is  that  the  method  employed  in  [8]  has  a  clear  cut  advantage. 


E.9  Case  IX 

In  this  section  the  results  for  a  case  similar  to  Case  VIII  are  presented. 
Everything  is  the  same  except  that  the  right  boundary  condition  is  held  fixed, 

Tn  *  0.1.  The  results  are  shown  in  Figures  IXa,  b,  and  c.  Until  about  t=l ,  the 
results  for  the  two  cases  are  practically  identical.  At  steady  state  the  results 
are  completely  different  as  seen  by  a  comparison  of  Figures  VIII  and  IX.  From 
about  t  =  1.1  to  t  =  1.3,  the  integration  of  Case  IX  proceeded  significantly  more 
efficiently.  In  Figure  IXa  the  T-profiles  for  t  =  1.2  and  t  *  1.5  are  practically 
indistinguishable  since  t  =  1.2  is  nearly  steady  state.  (The  +'s  correspond  to 
t  =  1.2  and  the  X's  correspond  to  t  =  1.5.)  One  interesting  aspect  of  the  final 
T-profile  is  that  the  minimum  mesh  spacing  is  0(0.003),  considerably  smaller  than 
the  minimum  of  0(0.007)  of  Case  VIII.  If  Case  IX  were  to  be  solved  by  the  methods 
of  [8],  the  advantage  of  [8]  over  the  present  method  would  not  be  so  clear  cut 
because  of  the  increased  number  of  mesh  points  required  to  obtain  a  mesh  size  of 
0.003  on  a  uniform  grid. 


1-59 


TEXTRON 


DT/DX 


0.40  0. 60  0.80 

T-SEC. 


1.00  1.20  1.40 


Figure  IXb  -  x-Trajectories  for  Case  IX 
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Figure  IXc  -  v-Trajectorles  and  Integration  Step  Sizes  for  Case  IX 

1-61 


TEXTRON 


Bell  Aerospace 


II.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  Conclusions 

During  this  study  a  variety  of  mesh  spacing  formulas  have  been  used.  The 
better  ones,  e.g..  Cases  VII,  VIII,  and  IX  did  a  reasonably  good  job  of  placing 
the  mesh  points  where  they  were  most  needed.  These  formulas  are  based  more  on 
experiment  and  intuition  than  on  theory.  Despite  this  lack  of  theoretical 
foundation,  the  variable  meshes  so  produced  are  far  better  than  uniform  meshes. 

For  example,  in  Case  VIII  31  mesh  points  produced  about  the  same  accuracy  as  200 
uniformly  spaced  points.  Of  course,  for  Burger's  equation  there  are  three 
dependent  variables  per  mesh  point  for  the  formulation  used  in  this  study,  versus 
one  per  mesh  point  for  the  usual  formulation,  e.g.,  reference  [8].  For  large 
partial  differential  systems,  e.g.,  those  involving  many  chemical  dependent 
variables,  the  ratio  (dependent  variables/mesh  points)  would  approach  2n,  where  n 
is  the  number  of  basic  dependent  variables.  Thus  if  a  variable  mesh  scheme  uses 
only  half  as  many  mesh  points  as  a  uniform  mesh  method  for  such  problems,  the 
former  method  would  be  competitive. 

The  foregoing  conclusion  is  true  providing  the  time  integration  can  be  performed 
as  efficiently  for  variable  meshes  as  it  is  for  uniform  meshes.  However,  this  was 
not  the  case.  In  Case  VIII  the  integration  statistics  were  much  worse  than  those 
for  a  similar  case  reported  in  [8],  One  experiment  performed  during  this  study  was 
to  take  a  variable  mesh  case  that  was  having  integration  difficulty  (e.g.,  small  and 
erratic  step  sizes,  build-up  of  spurious  oscillations,  difficulty  with  corrector 
convergence),  stop  the  integration  at  some  point  before  "hang  up,"  redefine  the 
mesh  spacing  formula  so  that  the  mesh  would  stay  fixed  from  then  on,  and  then  restart 
the  integration.  The  results  were  dramatic;  the  oscillations  died  out  and  the 
integration  proceeded  smoothly.  Of  course,  the  fixed  mesh  soon  became  inadequate 
to  describe  the  solution. 

What  is  the  reason  for  this  integration  difficulty?  Most  likely  it  is  due  to 
poor  conditioning,  especially  when  the  T-profiles  become  steep  and  have  regions  of 
high  curvature.  This  was  already  discussed  in  conjunction  with  Case  VIII.  Attempts 
were  made  to  solve  Case  VIII  with  8  =  0.003  (instead  of  3  =  0.005  that  was  eventually 
used).  Considerable  difficulty  in  obtaining  initial  condition  convergence  was 
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encountered.  Even  when  convergence  was  obtained,  the  distribution  was  peculiar 
and  |v.j  -  v.jexactl  *  was  relatively  large.  However,  the  residuals  of  the  algebraic 
equations  that  were  being  solved  had  been  driven  to  small  values.  This  phenomenon 
is  practically  the  definition  of  ill  conditioning. 


Another  indication  of  this  poor  conditioning  was  the  disappointingly  large 
ratios  (NRE/NSTEP)**  and  (NJE/NSTEP)**  that  occurred  in  most  of  our  numerical 
examples.  Such  large  ratios  are  symptomatic  of  ill -conditioning. 


B.  Recomnendations 

Although  definite  progress  has  been  made  in  developing  an  automatic  mesh  varying 
method,  certain  obstacles  {especially  the  conditioning  problem)  preclude  writing  a 
production  program  based  on  this  method.  Additional  analysis  and  experimentation  is 
required  in  the  following  areas: 

1)  Replace  the  banded  matrix  equation  solver  within  GEARIB  by  one  which  also 
provides  an  estimate  of  the  condition  number  of  the  matrix,  e.g.,  subroutine 
SGBCO  described  in  [9].  (This  change  may  require  extensive  coding  revisions 
of  GEARIB  to  accommodate  the  different  storage  arrangement  of  SGBCO.) 

Study  the  effect  of  the  mesh  spacing  formula  (including  any  "tuning"  constants) 
on  the  condition  of  the  matrix  used  for  the  quasi -Newton  corrector  iterations. 

2)  Continue  investigation  of  mesh  spacing  formulas.  Consider  introducing 
additional  mesh  point  dependent  variables  for  the  purpose  of  estimating 
third  and  fourth  derivatives  of  T  (with  respect  to  x)  in  order  to  compute 

a  mesh  spacing  formula  more  consistent  with  the  theoretical  error  estimates 
of  equation  (19). 


*Here  v..  exact  is  the  derivative  with  respect  to  x  of  equation  (39)  evaluated  at 
xi. 

**NSTEP:  Number  of  successful  integration  steps 

NRE:  Number  of  corrector  iterations 

NJE:  Number  of  Jacobian  evaluations 
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3)  From  the  standpoint  of  integration  efficiency  our  most  successful  example 
was  Case  V.  What  made  this  case  run  efficiently  was  that  the  mesh  points 
moved  with  the  wave  rather  than  passing  through  it.  Consequently,  once 
the  wave  established  itself,  the  dependent  variables  were  either  constant 
or  linear  in  time.  This  suggests  developing  methods  by  which  mesh  points 
are  added  to  (deleted  from)  the  field  as  they  are  needed  (unneeded).  By 
so  doing,  shifting  of  the  entire  set  of  mesh  points  just  to  transfer  a 
single  mesh  point  from  one  place  to  another  would  be  averted.  (This  was 
seen  most  dramatically  in  the  "crossover  effect"  of  Figures  1 1 b ,  1 1 1 b ,  and 
IVb,  but  it  is  also  evident  to  a  lesser  extent  in  Figures  lb,  Vllb,  VI 1 1 b , 
and  IXb.)  These  mesh  point  additions  (deletions)  might  be  done  at  discrete 
times,  in  which  case  the  method  resembles  the  regridding  approaches  of 
references  [5]  and  [6].  Preferably,  this  would  be  accomplished  in  a 
continuous  manner  by  allowing  a  single  mesh  point  to  split  in  two  or  by  two 
mesh  points  coalescing  into  one. 

4)  One  of  the  assumptions  throughout  this  study  was  that  x-j  would  coincide 
with  the  left  boundary  (say  xL)  of  the  mathematical  problem  and  x^  with 
the  right  boundary  (say  xR),  and  that  the  number  of  mesh  points  N  would 
be  constant.  According  to  formula  (10),  A  then  varies  as  the  solution 
evolves.  Now  the  spatial  truncation  error  is  related  to  A  (equation  (21)). 

A  preferable  approach  would  be  to  specify  the  desired  level  of  truncation 
error  (and  therefore  A)  and  then  determine  the  number  of  mesh  points  N  and/or 
the  length  of  the  finite  difference  mesh  [x-j ,  x^]  consistent  with  this 
specification.  If  x^  <  xL  <  X£  and/or  x^  -j  c  xR  <  x^,  the  boundary  conditions 
at  xL  and  xR  would  be  expressed  in  terms  of  values  at  two  consecutive  mesh 
points  instead  of  the  simpler  expressions  of  system  (13).  This  would  increase 
the  bandwidth  at  the  beginning  and  end  of  the  Jacobian  matrix.  However, 
special  formulas  could  be  developed  which  would  allow  solution  of  such 
modified  matrices  with  very  little  additional  cost.  A  simplification  arising 
out  of  such  a  formulation  is  that  A  is  not  a  variable,  and  so  the  modification 
of  the  GEARIB  banded  equation  solver  to  accommodate  a  possibly  full  last 
column  would  no  longer  be  required. 


1 1-3 


Bell  Aerospace 


5)  Ultimately,  it  is  desired  to  extend  these  variable  mesh  methods  to  gas 
dynamic  problems  involving  chemistry  and  radiation.  Even  "standard" 
gas  dynamic  problems  exhibit  phenomena  not  present  in  solutions  to 
Burger's  equation,  e.g.,  waves  traveling  at  variable  speeds,  waves  being 
reflected  from  boundaries,  waves  traveling  in  opposite  directions  and 
impinging  on  each  other.  Thus,  even  after  variable  mesh  methods  are 
perfected  for  Burger's  equation,  there  is  still  much  work  to  be  done 
before  these  methods  can  be  successfully  applied  to  realistic  problems. 
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